C************E02DAF.FOR************
      SUBROUTINE E02DAF(M, PX, PY, X, Y, F, W, LAMDA, MU, POINT,NPOINT,
     * DL, C, NC, WS, NWS, EPS, SIGMA, RANK, IFAIL)
C     MARK 6 RELEASE. NAG COPYRIGHT 1977
C
C     THIS IS A CALLING ROUTINE WHICH SETS UP WORK SPACE ARRAYS
C     FOR THE SUBROUTINE E02DAZ, AND THEN CALLS IT.
C
C     .. SCALAR ARGUMENTS ..
      REAL*4 EPS, SIGMA
      INTEGER IFAIL, M, NC, NPOINT, NWS, PX, PY, RANK
C     .. ARRAY ARGUMENTS ..
      REAL*4 C(NC), DL(NC), F(M), LAMDA(PX), MU(PY), W(M),
     * WS(NWS),X(M), Y(M)
      INTEGER POINT(NPOINT)
C     ..
C     .. LOCAL SCALARS ..
C$P 1
      DOUBLE PRECISION SRNAME
      INTEGER BW, IA, IB, ID, IDT, IL, IRT, IT, K
C     .. FUNCTION REFERENCES ..
      INTEGER P01AAF
C     .. SUBROUTINE REFERENCES ..
C     E02DAZ
C     ..
      DATA SRNAME /8H E02DAF /
      BW = 4 + 3*(PY-4)
      K = 0
      IF (NWS.LT.2*NC*(2+BW)+BW) K = 4
      IF (K.GT.0) GO TO 20
      ID = 1
      IB = ID + NC
      IT = IB + NC
      IA = IT + NC*BW
      IDT = IA + BW
      IL = IDT + NC
      IRT = IL + NC
      CALL E02DAZ(M, PX, PY, X, Y, F, W, LAMDA, MU, POINT, NPOINT,DL, C,
     * NC, WS(ID), WS(IB), WS(IT), WS(IA), BW, WS(IDT),WS(IL), WS(IRT),
     * EPS, SIGMA, RANK, K)
      IF (K.NE.0) GO TO 20
      IFAIL = K
      RETURN
   20 IFAIL = P01AAF(IFAIL,K,SRNAME)
      RETURN
      END
      SUBROUTINE E02DAZ(M, PX, PY, X, Y, F, W, LAMDA, MU, POINT,NPOINT,
     * DL, C, NC, DIAG, B, T, A, BW, DTILDE, LAM, RTILDE,EPS, SIGMA,
     * RANK, IFAIL)
C     MARK 6 RELEASE. NAG COPYRIGHT 1977
C     MARK 7 REVISED IER-142 (DEC 1978)
C     MARK 10A REVISED. IER-395 (OCT 1982).
C
C     THIS SUBROUTINE FORMS THE MINIMAL, WEIGHTED LEAST-SQUARES
C     BICUBIC SPLINE SURFACE FIT (HAYES, J.G. AND HALLIDAY,J., THE
C     LEAST-SQUARES FITTING OF CUBIC SPLINE SURFACES TO GENERAL
C     DATA SETS, NPL DIVISIONAL REPORT DNAC 22 (1972)) WITH
C     PRESCRIBED NON-DECREASING KNOTS LAMDA(Q), Q= 1(1)PX AND
C     MU(R), R= 1(1)PY TO THE M DATA POINTS X(L), Y(L), F(L),
C     L= 1(1)M.   THE FORM OF THE FIT IS
C     F(X,Y)= C(I,J)*MI(X)*NJ(Y),
C     SUMMED FOR I= 1,2,...PX-4 AND J= 1,2,...PY-4, WHERE MI(X) AND
C     NJ(Y) ARE THE NORMALISED B-SPLINES DEFINED ON THE KNOT SETS
C     LAMDA AND MU RESPECTIVELY. (DEBOOR, C., ON CALCULATING WITH
C     B-SPLINES, J. APPROX. THEORY, VOL 6 (1972), COX, M.G., THE
C     NUMERICAL EVALUATION OF B-SPLINES, NPL DIVISIONAL REPORT DNAC
C     4 (1971)). THE KNOTS WITH Q= 5,6,...,PX-4 AND R= 5,6,...,PY-4
C     SHOULD BE CHOSEN INTERIOR TO THE X AND Y DATA VALUES, RESP.
C     THE REMAINING KNOTS, FOUR AT EACH END OF EACH DATA RANGE, ARE
C     CHOSEN BY THE SUBROUTINE TO COINCIDE AT THE APPROPRIATE
C     EXTREME DATA VALUE OF X OR Y. TWO, OR THREE, COINCIDENT KNOTS
C     INTERIOR TO THE RANGE OF AN INDEPENDENT VARIABLE CAUSE LOSS
C     OF CONTINUITY OF THE FIT IN, RESPECTIVELY, THE SECOND AND
C     FIRST DERIVATIVE WITH RESPECT TO THAT VARIABLE. THE ARRAY
C     POINT CONTAINS INDEXING INFORMATION SO THAT THE DATA POINTS
C     CAN BE ACCESSED IN SUCH AN ORDER THAT, IF THE (X,Y)-PLANE IS
C     THOUGHT OF AS BEING DIVIDED INTO RECTANGULAR PANELS BY THE
C     TWO SETS OF KNOTS, ALL DATA IN A PANEL OCCURS BEFORE DATA IN
C     SUCCEEDING PANELS, WHERE THE PANELS ARE NUMBERED FROM BOTTOM
C     TO TOP AND THEN LEFT TO RIGHT WITH THE USUAL ARRANGEMENT OF
C     AXES. A DATA POINT AT A KNOT IS CONSIDERED TO BE IN THE
C     HIGHER NUMBERED PANEL, WITH OBVIOUS EXCEPTIONS AT THE UPPER
C     ENDS OF THE RANGES OF THE INDEPENDENT VARIABLES. A FORTRAN
C     SUBROUTINE E02ZAF IS AVAILABLE TO SET UP POINT IN THIS WAY.
C     THE OBSERVATION MATRIX IS REDUCED TO UPPER TRIANGULAR FORM BY
C     GIVENS TRANSFORMATIONS WITHOUT SQUARE ROOTS (GENTLEMAN,
C     W.M., LEAST SQUARES COMPUTATIONS BY GIVENS TRANSFORMATIONS
C     WITHOUT SQUARE ROOTS, J. INST. MATH. APPLICS. VOL 12 NO 3
C     (1973)), EACH ROW BEING ROTATED INTO THE TRIANGLE AS IT IS
C     FORMED, THUS AVOIDING THE NEED TO STORE THE ENTIRE
C     OBSERVATION MATRIX. THE RANK OF THE SYSTEM IS DETERMINED AS
C     THE NUMBER OF DIAGONAL ELEMENTS OF THE REDUCED TRIANGLE WHOSE
C     SQUARES, DIVIDED BY THE MEAN SQUARED WEIGHT, ARE NOT LESS
C     THAN EPS. IF THE NUMBER OF DECIMAL DIGITS IN THE COMPUTER
C     REPRESENTATION OF A REAL NUMBER IS T, THEN 10**(-T) IS A
C     SUITABLE VALUE FOR EPS IN MOST PRACTICAL APPLICATIONS.
C     THE SUBROUTINE CHECKS THAT
C     1) EACH SET OF KNOTS IS IN NON-DECREASING ORDER,
C     2) NO KNOT HAS MULTIPLICITY GREATER THAN 4,
C     3) THE DATA IS ORDERED IN PANELS AS DESCRIBED ABOVE,
C     4) MY>1, PX AND PY ARE EACH EQUAL TO AT LEAST 8,
C     NC = (PX-4)*(PY-4) AND NWS AND NPOINT ARE LARGE ENOUGH.
C     5) NOT ALL THE WEIGHTS ARE ZERO
C     IF ANY OF THE CHECKS FAILS THE SUBROUTINE EXITS IMMEDIATELY
C     WITH THE INTEGER IFAIL SET EQUAL TO THE NUMBER OF THE CHECK
C     WHICH HAS FAILED.   ON NORMAL EXIT IFAIL EQUALS ZERO.
C
C     .. SCALAR ARGUMENTS ..
      REAL*4 EPS, SIGMA
      INTEGER BW, IFAIL, M, NC, NPOINT, PX, PY, RANK
C     .. ARRAY ARGUMENTS ..
      REAL*4 A(BW), B(NC), C(NC), DIAG(NC), DL(NC), DTILDE(NC)
     *, F(M),LAM(NC), LAMDA(PX), MU(PY), RTILDE(NC,BW), T(NC,BW), W(M),
     *X(M), Y(M)
      INTEGER POINT(NPOINT)
C     ..
C     .. LOCAL SCALARS ..
      REAL*4 BCOL, BROW, COSINE, D, DPRIME, DX4, DX5, DX6,
     * DX7, DX8,DX9, DY4, DY5, DY6, DY7, DY8, DY9, E2, E3, E4, E5,
     * LAMDA1,LAMDA2, LAMDA3, LAMDA4, LAMDA5, LAMDA6, M1, M2, M3, MNW,
     *MU1, MU2, MU3, MU4, MU5, MU6, MXL, ROWEL, S, SINE, SQ, TCOL,TROW,
     * WI, XI, YI
      INTEGER BW1, BWMINL, I, IADRES, II, IITEMP, ITEMP, J, JJ,JPLUSL,
     * JX, JXOLD, JY, K4, K7, K, L4, L7, L, LL, LPLUSU, NP,R, ROWNUM,
     * S1, S2, T2, U, WN
C     .. LOCAL ARRAYS ..
      REAL*4 MX(4), MY(4)
C     ..
      K4 = PX - 4
      L4 = PY - 4
      K7 = PX - 7
      L7 = PY - 7
      NP = K7*L7
      BW1 = BW - 1
      IFAIL = 4
      IF (M.LE.1 .OR. PX.LT.8 .OR. PY.LT.8 .OR. NC.NE.K4*L4 .OR.NPOINT.
     *LT.M+NP) RETURN
C
C     SET THE FOUR KNOTS AT  EACH END OF EACH DATA
C     RANGE EQUAL TO THE APPROPRIATE EXTREME DATA VALUE.
C
      M1 = X(1)
      M2 = M1
      DO 20 I=1,M
         XI = X(I)
         IF (XI.GT.M1) M1 = XI
         IF (XI.LT.M2) M2 = XI
   20 CONTINUE
      DO 40 I=1,4
         LAMDA(I) = M2
         ITEMP = PX - I + 1
         LAMDA(ITEMP) = M1
   40 CONTINUE
      M1 = Y(1)
      M2 = M1
      DO 60 I=1,M
         YI = Y(I)
         IF (YI.GT.M1) M1 = YI
         IF (YI.LT.M2) M2 = YI
   60 CONTINUE
      DO 80 I=1,4
         MU(I) = M2
         ITEMP = PY - I + 1
         MU(ITEMP) = M1
   80 CONTINUE
C
C     CHECK THAT EACH KNOT SET IS IN NON-DECREASING ORDER.
C
      IFAIL = 1
      DO 100 I=2,PX
         IF (LAMDA(I-1).GT.LAMDA(I)) RETURN
  100 CONTINUE
      DO 120 I=2,PY
         IF (MU(I-1).GT.MU(I)) RETURN
  120 CONTINUE
C
C     CHECK THAT NO KNOT HAS MULTIPLICITY GREATER THAN FOUR
C
      IFAIL = 2
      DO 140 I=1,K4
         IF (LAMDA(I).EQ.LAMDA(I+4)) RETURN
  140 CONTINUE
      DO 160 I=1,L4
         IF (MU(I).EQ.MU(I+4)) RETURN
  160 CONTINUE
C
C     CHECK THAT NOT ALL THE WEIGHTS ARE ZERO
C
      IFAIL = 5
      DO 180 I=1,M
         IF (W(I).NE.0.0) IFAIL = 0
  180 CONTINUE
      IF (IFAIL.NE.0) RETURN
      IFAIL = 3
C
C     INITIALISE UPPER TRIANGLE AND DIAGONAL TO ZERO
C
      DO 220 I=1,NC
         DIAG(I) = 0
         B(I) = 0
         DO 200 J=2,BW
            T(I,J) = 0
  200    CONTINUE
  220 CONTINUE
      SQ = 0
      JXOLD = 0
      MNW = 0
      S1 = 0
C
C     ROWNUM RECORDS THE COLUMN NUMBER OF THE FIRST
C     NON-ZERO ELEMENT IN THIS ROW OF THE OBSERVATION MATRIX,
C     I.E. THE NUMBER OF THE FIRST ROW OF THE TRIANGLE WITH
C     WHICH IT IS TO BE ROTATED
C
      DO 440 I=1,NP
         IADRES = I + M
         IF (POINT(IADRES).LE.0) GO TO 440
         JX = (I-1)/L7
         JY = I - JX*L7
         JX = JX + 1
         ROWNUM = (JX-1)*L4 + JY
         IF (JX.EQ.JXOLD) GO TO 240
         LAMDA1 = LAMDA(JX+1)
         LAMDA2 = LAMDA(JX+2)
         LAMDA3 = LAMDA(JX+3)
         LAMDA4 = LAMDA(JX+4)
         LAMDA5 = LAMDA(JX+5)
         LAMDA6 = LAMDA(JX+6)
         DX4 = 1.0/(LAMDA4-LAMDA1)
         DX5 = 1.0/(LAMDA5-LAMDA2)
         DX6 = 1.0/(LAMDA6-LAMDA3)
         DX7 = 1.0/(LAMDA4-LAMDA2)
         DX8 = 1.0/(LAMDA5-LAMDA3)
         DX9 = 1.0/(LAMDA4-LAMDA3)
         JXOLD = JX
  240    MU1 = MU(JY+1)
         MU2 = MU(JY+2)
         MU3 = MU(JY+3)
         MU4 = MU(JY+4)
         MU5 = MU(JY+5)
         MU6 = MU(JY+6)
         DY4 = 1.0/(MU4-MU1)
         DY5 = 1.0/(MU5-MU2)
         DY6 = 1.0/(MU6-MU3)
         DY7 = 1.0/(MU4-MU2)
         DY8 = 1.0/(MU5-MU3)
         DY9 = 1.0/(MU4-MU3)
  260    IADRES = POINT(IADRES)
         IF (IADRES.LE.0) GO TO 440
C
C     INITIALISE CURRENT ROW OF OBSERVATION MATRIX TO ZERO
C
         DO 280 J=1,BW
            A(J) = 0
  280    CONTINUE
C
C     FETCH CURRENT DATA POINT XI, YI
C
         WI = W(IADRES)**2
         IF (WI.EQ.0.00) GO TO 260
         XI = X(IADRES)
         YI = Y(IADRES)
         BROW = F(IADRES)
         S1 = S1 + 1
         MNW = MNW + WI
C
C     CHECK THAT DATA POINTS ARE ACCESSED IN PANEL ORDER
C
         IF (XI.LT.LAMDA3 .OR. XI.GT.LAMDA4 .OR. (XI.EQ.LAMDA4.AND. JX.
     *   NE.K7) .OR. YI.LT.MU3 .OR. YI.GT.MU4 .OR.(YI.EQ.MU4 .AND. JY.
     *   NE.L7)) RETURN
C
C     EVALUATE THE FOUR NON-ZERO CUBIC B-SPLINES AT XI,
C     AND STORE IN MX(1) TO MX(4)
C
         E5 = LAMDA5 - XI
         E4 = LAMDA4 - XI
         E3 = XI - LAMDA3
         E2 = XI - LAMDA2
         M2 = E3*DX8*DX9
         M1 = E4*DX7*DX9
         M3 = E3*M2*DX6
         M2 = (E2*M1+E5*M2)*DX5
         M1 = E4*M1*DX4
         MX(4) = E3*M3
         MX(3) = E2*M2 + (LAMDA6-XI)*M3
         MX(2) = (XI-LAMDA1)*M1 + E5*M2
         MX(1) = E4*M1
C
C     EVALUATE THE FOUR NON-ZERO CUBIC B-SPLINES AT YI,
C     AND STORE IN MY(1) TO MY(4)
C
         E5 = MU5 - YI
         E4 = MU4 - YI
         E3 = YI - MU3
         E2 = YI - MU2
         M2 = E3*DY8*DY9
         M1 = E4*DY7*DY9
         M3 = E3*M2*DY6
         M2 = (E2*M1+E5*M2)*DY5
         M1 = E4*M1*DY4
         MY(4) = E3*M3
         MY(3) = E2*M2 + (MU6-YI)*M3
         MY(2) = (YI-MU1)*M1 + E5*M2
         MY(1) = E4*M1
         K = 0
C
C     EVALUATE ELEMENTS OF THIS ROW OF OBSERVATION MATRIX
C     WHICH ARE BICUBIC B-SPLINES FORMED AS PRODUCTS OF
C     CUBIC B-SPLINES IN XI AND YI RESPECTIVELY
C
         DO 320 L=1,4
            MXL = MX(L)
            DO 300 J=1,4
               IITEMP = K + J
               A(IITEMP) = MXL*MY(J)
  300       CONTINUE
            K = K + L4
  320    CONTINUE
C
C     ROTATE THIS ROW OF OBSERVATION MATRIX INTO THE
C     TRIANGLE BY GIVENS TRANSFORMATIONS WITHOUT SQUARE ROOTS
C
         J = ROWNUM
C
C     L INDEXES ELEMENTS IN BAND
C
         DO 420 LL=1,BW
            L = LL - 1
C
C     NO ROTATION IS PERFORMED IF CURRENT ELEMENT
C     OF OBSERVATION MATRIX IS ZERO
C
            IF (WI.EQ.0.00) GO TO 260
            ROWEL = A(L+1)
            S = WI*ROWEL
            IF (S*ROWEL.EQ.0.0) GO TO 420
            JPLUSL = J + L
            BWMINL = BW - L
            D = DIAG(JPLUSL)
            IF (D.NE.0.00) GO TO 340
            DIAG(JPLUSL) = S*ROWEL
            COSINE = 0
            SINE = 1.0/ROWEL
            GO TO 360
  340       DPRIME = D + S*ROWEL
            DIAG(JPLUSL) = DPRIME
            COSINE = D/DPRIME
            SINE = S/DPRIME
  360       WI = WI*COSINE
            IF (BWMINL.LT.2) GO TO 400
            DO 380 U=2,BWMINL
               LPLUSU = L + U
               TCOL = T(JPLUSL,U)
               TROW = A(LPLUSU)
               T(JPLUSL,U) = COSINE*TCOL + SINE*TROW
               A(LPLUSU) = TROW - ROWEL*TCOL
  380       CONTINUE
C
C     APPLY SAME TRANSFORMATIONS TO RIGHT HAND SIDE
C
  400       BCOL = B(JPLUSL)
            B(JPLUSL) = COSINE*BCOL + SINE*BROW
            BROW = BROW - ROWEL*BCOL
  420    CONTINUE
C
C     ADD CONTRIBUTION OF THIS ROW TO SUM OF SQUARES
C     OF RESIDUAL RIGHT HAND SIDES
C
         SQ = SQ + WI*BROW**2
         GO TO 260
  440 CONTINUE
C
C     FORM MEAN SQUARED WEIGHT
C
      M1 = S1
      MNW = MNW/M1
      R = 0
C
C     DETERMINE RANK DEFICIENCY R AS NUMBER OF
C     SUFFICIENTLY SMALL DIAGONAL ELEMENTS OF REDUCED
C     TRIANGLE, AND ROTATE THE CORRESPONDING ROW, WITH
C     THE DIAGONAL TERM MADE ZERO, INTO THE REMAINDER
C     OF THE TRIANGLE
C
      DO 660 I=1,NC
         WI = DIAG(I)
         DL(I) = WI/MNW
         IF (WI.GT.EPS*MNW) GO TO 660
         DIAG(I) = 0
         R = R + 1
         BROW = B(I)
         DO 460 L=2,BW
            A(L-1) = T(I,L)
  460    CONTINUE
         A(BW) = 0
         IF (I.EQ.NC) GO TO 640
         II = I + 1
         DO 620 L=II,NC
            IF (WI.EQ.0.00) GO TO 660
            ROWEL = A(1)
            S = WI*ROWEL
            WN = NC - L + 1
            IF (WN.GT.BW) WN = BW
            IF (S*ROWEL.EQ.0.00) GO TO 560
            D = DIAG(L)
            IF (D.NE.0.00) GO TO 480
            DIAG(L) = S*ROWEL
            COSINE = 0
            SINE = 1.0/ROWEL
            GO TO 500
  480       DPRIME = D + S*ROWEL
            DIAG(L) = DPRIME
            COSINE = D/DPRIME
            SINE = S/DPRIME
  500       WI = WI*COSINE
            IF (WN.LT.2) GO TO 540
            DO 520 U=2,WN
               TCOL = T(L,U)
               TROW = A(U)
               T(L,U) = COSINE*TCOL + SINE*TROW
               A(U-1) = TROW - ROWEL*TCOL
  520       CONTINUE
  540       BCOL = B(L)
            B(L) = COSINE*BCOL + SINE*BROW
            BROW = BROW - ROWEL*BCOL
            GO TO 600
  560       DO 580 U=2,WN
               A(U-1) = A(U)
  580       CONTINUE
  600       A(WN) = 0
  620    CONTINUE
  640    SQ = SQ + WI*BROW**2
  660 CONTINUE
      RANK = NC - R
C
C     IF MATRIX IS RANK DEFICIENT, FIND MINIMUM NORM SOLUTION
C
      IFAIL = 5
      IF (RANK.EQ.0) RETURN
      IF (R.EQ.0) GO TO 1240
C
C     FORM IN RTILDE THE TRANSPOSE OF THE MATRIX
C     OBTAINED FROM THE REDUCED TRIANGLE BY REMOVING ROWS AND
C     COLUMNS WITH SUFFICIENTLY SMALL DIAGONAL ELEMENTS, AND
C     THEN REVERSING THE ORDERING OF THE REMAINING ROWS AND
C     COLUMNS
C
      DO 700 I=1,RANK
         DTILDE(I) = 1.0
         DO 680 J=2,BW
            RTILDE(I,J) = 0
  680    CONTINUE
  700 CONTINUE
      S1 = 0
      WN = 0
      ITEMP = NC - 1
      DO 740 II=1,ITEMP
         I = NC - II
         IF (DIAG(I).EQ.0.00) GO TO 740
         WN = II + 1
         IF (WN.GT.BW) WN = BW
         S2 = S1 + 1
         S1 = S2
         T2 = 2
         DO 720 J=2,WN
            IITEMP = I + J - 1
            IF (DIAG(IITEMP).EQ.0.00) GO TO 720
            RTILDE(S2,T2) = T(I,J)
            S2 = S2 - 1
            T2 = T2 + 1
  720    CONTINUE
         IF (S2.EQ.S1) S1 = S1 - 1
  740 CONTINUE
C
C     FORM NEW RIGHT HAND SIDE IN LAM BY REMOVING
C     ELEMENTS OF OLD RIGHT HAND SIDE WHICH CORRESPOND TO A
C     DELETED ROW AND REVERSING ORDERING OF REMAINING ELEMENTS
C
      S1 = 0
      DO 760 II=1,NC
         I = NC + 1 - II
         IF (DIAG(I).EQ.0.00) GO TO 760
         S1 = S1 + 1
         LAM(S1) = B(I)
  760 CONTINUE
C
C     FORM SUCCESSIVELY IN A THOSE ROWS OF TRANSPOSE
C     OF REDUCED TRIANGLE WHICH HAD SUFFICIENTLY SMALL
C     DIAGONAL ELEMENTS
C
      WN = 0
      DO 770 I =1,NC
         IF (DIAG(I).EQ.0.00) GO TO 770
         II = I
         GO TO 775
  770 CONTINUE
  775 S1 = RANK + II
      DO 960 I=II,ITEMP
         IF (DIAG(I+1).NE.0.00) GO TO 960
         S2 = 1
         WN = BW
         IF (I.LT.BW) WN = I + 1
         DO 780 J=2,WN
            IITEMP = I - J + 2
            IF (DIAG(IITEMP).EQ.0.00) GO TO 780
            A(S2) = T(IITEMP,J)
            S2 = S2 + 1
  780    CONTINUE
         DO 800 J=S2,BW
            A(J) = 0
  800    CONTINUE
C
C     ROTATE THIS ROW INTO RTILDE BY GIVENS
C     TRANSFORMATIONS WITHOUT SQUARE ROOTS
C
         J = S1 - I
         S1 = S1 + 1
         WI = 1.0
         DO 940 L=J,RANK
            IF (WI.EQ.0.0) GO TO 960
            ROWEL = A(1)
            S = WI*ROWEL
            IF (S*ROWEL.EQ.0.00) GO TO 880
            WN = RANK - L + 1
            IF (WN.GT.BW) WN = BW
            D = DTILDE(L)
            IF (D.NE.0.00) GO TO 820
            DTILDE(L) = S*ROWEL
            COSINE = 0
            SINE = 1.0/ROWEL
            GO TO 840
  820       DPRIME = D + S*ROWEL
            DTILDE(L) = DPRIME
            COSINE = D/DPRIME
            SINE = S/DPRIME
  840       WI = WI*COSINE
            IF (WN.LT.2) GO TO 920
            DO 860 U=2,WN
               TCOL = RTILDE(L,U)
               TROW = A(U)
               RTILDE(L,U) = COSINE*TCOL + SINE*TROW
               A(U-1) = TROW - ROWEL*TCOL
  860       CONTINUE
            GO TO 920
  880       WN = RANK - L + 1
            IF (WN.GT.BW) WN = BW
            IF (WN.LT.2) GO TO 920
            DO 900 U=2,WN
               A(U-1) = A(U)
  900       CONTINUE
  920       A(WN) = 0
  940    CONTINUE
  960 CONTINUE
      L = 0
C
C     FORWARD SUBSTITUTION
C
      DO 1020 I=1,RANK
         IF (L.LT.BW) L = I
         S = LAM(I)
         S1 = I
         IF (L.LT.2) GO TO 1000
         DO 980 J=2,L
            S1 = S1 - 1
            S = S - RTILDE(S1,J)*LAM(S1)
  980    CONTINUE
 1000    LAM(I) = S
 1020 CONTINUE
C
C     BACKWARD SUBSTITUTION
C
      L = 0
      DO 1080 JJ=1,RANK
         J = RANK + 1 - JJ
         IF (L.LT.BW) L = L + 1
         S1 = J - 1
         S = LAM(J)/DTILDE(J)
         IF (L.LT.2) GO TO 1060
         DO 1040 I=2,L
            IITEMP = I + S1
            S = S - RTILDE(J,I)*LAM(IITEMP)
 1040    CONTINUE
 1060    LAM(J) = S
 1080 CONTINUE
C
C     PREMULTIPLY LAM BY RTILDE
C
      S1 = RANK + 1
      WN = 0
      DO 1160 I=1,NC
         S2 = S1
         S = 0
         IF (DIAG(I).EQ.0.0) GO TO 1100
         S1 = S1 - 1
         S = LAM(S1)
 1100    IF (WN.LT.BW) WN = I
         IF (WN.LT.2) GO TO 1140
         DO 1120 J=2,WN
            IITEMP = I - J + 1
            IF (DIAG(IITEMP).EQ.0.00) GO TO 1120
            S = S + T(IITEMP,J)*LAM(S2)
            S2 = S2 + 1
 1120    CONTINUE
 1140    C(I) = S
 1160 CONTINUE
C
C     ADD TO RESIDUAL SUM OF SQUARES THE CONTRIBUTION
C     OF DISCARDED SMALL DIAGONAL ELEMENTS OF THE REDUCED
C     TRIANGLE
C
      DO 1220 I=1,NC
         IF (DIAG(I).NE.0.00) GO TO 1220
         S = B(I)
         IF (I.EQ.NC) GO TO 1200
         WN = NC - I + 1
         IF (WN.GT.BW) WN = BW
         DO 1180 J=2,WN
            ITEMP = I + J - 1
            S = S - C(ITEMP)*T(I,J)
 1180    CONTINUE
 1200    SQ = SQ + MNW*DL(I)*C(I)*(C(I)-2.00*S)
 1220 CONTINUE
      GO TO 1320
C
C     BACKWARD SUBSTITUTION FOR FULL RANK CASE
C
 1240 L = 0
      DO 1300 JJ=1,NC
         J = NC + 1 - JJ
         IF (L.LT.BW) L = L + 1
         U = J - 1
         S = B(J)
         IF (L.LT.2) GO TO 1280
         DO 1260 I=2,L
            IITEMP = I + U
            S = S - T(J,I)*C(IITEMP)
 1260    CONTINUE
 1280    C(J) = S
 1300 CONTINUE
 1320 IFAIL = 0
      SIGMA = SQ
      RETURN
      END
C************E02DBF.FOR************
      SUBROUTINE E02DBF(M, PX, PY, X, Y, FF, LAMDA, MU, POINT,NPOINT, C,
     * NC, IFAIL)
C     MARK 6 RELEASE. NAG COPYRIGHT 1977
C
C     THIS SUBROUTINE EVALUATES AT THE POINTS WHOSE
C     CO-ORDINATES ARE GIVEN IN X AND Y A BICUBIC SPLINE
C     REPRESENTED
C     IN THE FORM PRODUCED BY E02DAF,  THE ARRAYS LAMBDA
C     AND MU SPECIFY THE COMPLETE SET OF KNOTS ASSOCIATED
C     WITH THE VARIABLES X AND Y, RESPECTIVELY.
C
C     THE SUBROUTINE CHECKS THAT
C     1) EACH SET OF KNOTS IS IN NON-DECREASING ORDER,
C     2) NO KNOT HAS MULTIPLICITY GREATER THAN 4,
C     3) THE DATA IS ORDERED IN PANELS AS DESCRIBED ABOVE,
C     4) M>0 AND PX AND PY ARE EACH EQUAL TO AT LEAST 8,
C     NC = (PX-4)*(PY-4) AND NPOINT IS LARGE ENOUGH,
C     5) NO POINT (X, Y) LIES OUTSIDE THE RECTANGLE
C     DEFINED BY THE EXTREME KNOT VALUES.
C     IF ANY OF THE CHECKS FAILS, THE SUBROUTINE EXITS IMMEDIATELY
C     WITH THE INTEGER IFAIL SET EQUAL TO THE NUMBER OF THE CHECK
C     WHICH HAS FAILED.   ON NORMAL EXIT IFAIL EQUALS ZERO.
C
C     .. SCALAR ARGUMENTS ..
      INTEGER IFAIL, M, NC, NPOINT, PX, PY
C     .. ARRAY ARGUMENTS ..
      REAL*4 C(NC), FF(M), LAMDA(PX), MU(PY), X(M), Y(M)
      INTEGER POINT(NPOINT)
C     ..
C     .. LOCAL SCALARS ..
C$P 1
      DOUBLE PRECISION SRNAME
      REAL*4 DX4, DX5, DX6, DX7, DX8, DX9, DY4, DY5, DY6, DY7,
     * DY8,DY9, E2, E3, E4, E5, LAMDA1, LAMDA2, LAMDA3, LAMDA4, LAMDA5,
     *LAMDA6, M1, M2, M3, MU1, MU2, MU3, MU4, MU5, MU6, MXL, TEMP,WI,
     * XI, YI
      INTEGER I, IADRES, ITEMP, JJ, JX, JXOLD, JY, K4, K7, K, KK,L4, L7,
     * LL, NP
C     .. LOCAL ARRAYS ..
      REAL*4 MX(4), MY(4)
C     .. FUNCTION REFERENCES ..
      INTEGER P01AAF
C     ..
      DATA SRNAME /8H E02DBF /
      K = 0
      K7 = PX - 7
      L7 = PY - 7
      NP = K7*L7
      K4 = PX - 4
      L4 = PY - 4
      IF (M.LE.0 .OR. PX.LT.8 .OR. PY.LT.8 .OR. NC.NE.K4*L4 .OR.NPOINT.
     *LT.M+NP) K = 4
      IF (K.NE.0) GO TO 220
C
C     CHECK THAT EACH KNOT SET IS IN NON-DECREASING ORDER.
C
      DO 20 I=2,PX
         IF (LAMDA(I-1).GT.LAMDA(I)) K = 1
   20 CONTINUE
      DO 40 I=2,PY
         IF (MU(I-1).GT.MU(I)) K = 1
   40 CONTINUE
      IF (K.NE.0) GO TO 220
C
C     CHECK THAT NO KNOT HAS MULTIPLICITY GREATER THAN FOUR
C
      DO 60 I=1,K4
         IF (LAMDA(I).EQ.LAMDA(I+4)) K = 2
   60 CONTINUE
      DO 80 I=1,L4
         IF (MU(I).EQ.MU(I+4)) K = 2
   80 CONTINUE
      IF (K.NE.0) GO TO 220
C
C     CHECK THAT DATA POINTS ARE NOT OUTSIDE THE RECTANGLE
C     DEFINED BY THE EXTREME KNOTS.
C
      LAMDA1 = LAMDA(4)
      LAMDA2 = LAMDA(PX-3)
      MU1 = MU(4)
      MU2 = MU(PY-3)
      DO 100 I=1,M
         XI = X(I)
         YI = Y(I)
         IF (XI.LT.LAMDA1 .OR. XI.GT.LAMDA2 .OR. YI.LT.MU1 .OR.YI.GT.
     *   MU2) K = 5
         IF (K.NE.0) GO TO 220
  100 CONTINUE
      JXOLD = 0
      DO 200 I=1,NP
         IADRES = I + M
         IF (POINT(IADRES).LE.0) GO TO 200
         JX = (I-1)/L7
         JY = I - JX*L7
         JX = JX + 1
         IF (JX.EQ.JXOLD) GO TO 120
         LAMDA1 = LAMDA(JX+1)
         LAMDA2 = LAMDA(JX+2)
         LAMDA3 = LAMDA(JX+3)
         LAMDA4 = LAMDA(JX+4)
         LAMDA5 = LAMDA(JX+5)
         LAMDA6 = LAMDA(JX+6)
         DX4 = 1.0/(LAMDA4-LAMDA1)
         DX5 = 1.0/(LAMDA5-LAMDA2)
         DX6 = 1.0/(LAMDA6-LAMDA3)
         DX7 = 1.0/(LAMDA4-LAMDA2)
         DX8 = 1.0/(LAMDA5-LAMDA3)
         DX9 = 1.0/(LAMDA4-LAMDA3)
         JXOLD = JX
  120    MU1 = MU(JY+1)
         MU2 = MU(JY+2)
         MU3 = MU(JY+3)
         MU4 = MU(JY+4)
         MU5 = MU(JY+5)
         MU6 = MU(JY+6)
         DY4 = 1.0/(MU4-MU1)
         DY5 = 1.0/(MU5-MU2)
         DY6 = 1.0/(MU6-MU3)
         DY7 = 1.0/(MU4-MU2)
         DY8 = 1.0/(MU5-MU3)
         DY9 = 1.0/(MU4-MU3)
  140    IADRES = POINT(IADRES)
         IF (IADRES.LE.0) GO TO 200
         XI = X(IADRES)
         YI = Y(IADRES)
C
C     CHECK THAT DATA POINTS ARE ACCESSED IN PANEL ORDER
C
         IF (XI.LT.LAMDA3 .OR. XI.GT.LAMDA4 .OR. (XI.EQ.LAMDA4.AND. JX.
     *   NE.K7) .OR. YI.LT.MU3 .OR. YI.GT.MU4 .OR.(YI.EQ.MU4 .AND. JY.
     *   NE.L7)) K = 3
         IF (K.NE.0) GO TO 220
C
C     EVALUATE THE FOUR NON-ZERO CUBIC B-SPLINESAT XI,
C     AND STORE IN MX(1) TO MX(4).
C
         E5 = LAMDA5 - XI
         E4 = LAMDA4 - XI
         E3 = XI - LAMDA3
         E2 = XI - LAMDA2
         M2 = E3*DX8
         M1 = E4*DX7
         M3 = E3*M2*DX6
         M2 = (E2*M1+E5*M2)*DX5
         M1 = E4*M1*DX4
         WI = DX9
         MX(4) = WI*E3*M3
         MX(3) = WI*(E2*M2+(LAMDA6-XI)*M3)
         MX(2) = WI*((XI-LAMDA1)*M1+E5*M2)
         MX(1) = WI*E4*M1
C
C     EVALUATE THE FOUR NON-ZERO CUBIC B-SPLINES AT YI,
C     AND STORE IN MY(1) TO MY(4).
C
         E5 = MU5 - YI
         E4 = MU4 - YI
         E3 = YI - MU3
         E2 = YI - MU2
         M2 = E3*DY8
         M1 = E4*DY7
         M3 = E3*M2*DY6
         M2 = (E2*M1+E5*M2)*DY5
         M1 = E4*M1*DY4
         WI = DY9
         MY(4) = WI*E3*M3
         MY(3) = WI*(E2*M2+(MU6-YI)*M3)
         MY(2) = WI*((YI-MU1)*M1+E5*M2)
         MY(1) = WI*E4*M1
C
C     EVALUATE THE BICUBIC SPLINE FUNCTION.
C
         KK = (JX-1)*L4 + JY
         TEMP = 0
         DO 180 LL=1,4
            MXL = MX(LL)
            DO 160 JJ=1,4
               ITEMP = KK + JJ - 1
               TEMP = TEMP + MXL*MY(JJ)*C(ITEMP)
  160       CONTINUE
            KK = KK + L4
  180    CONTINUE
         FF(IADRES) = TEMP
         GO TO 140
  200 CONTINUE
      IFAIL = K
      RETURN
  220 IFAIL = P01AAF(IFAIL,K,SRNAME)
      RETURN
      END

C************E02ZAF.FOR************
      SUBROUTINE E02ZAF(PX, PY, LAMDA, MU, M, X, Y, POINT, NPOINT,ADRES,
     * NADRES, IFAIL)
C     MARK 6 RELEASE. NAG COPYRIGHT 1977
C     MARK 10 REVISED. IER-385 (JUN 1982).
C
C     THE REAL ARRAYS LAMDA AND MU OF DIMENSION AT LEAST
C     PX-4 AND PY-4 RESPECTIVELY EACH CONTAIN A SEQUENCE
C     OF NON-DECREASING VALUES.  THE POINTS WITH COORDINATES
C     GIVEN IN X(1) TO X(M) AND Y(1) TO Y(M) ARE SORTED INTO
C     THE RECTANGULAR PANELS DETERMINED BY THE LAMDA AND MU
C     VALUES.  ARRAYS X AND Y ARE UNCHANGED, THE ORDERING
C     BEING INDICATED BY NADRES LINKED LISTS STORED IN ARRAY
C     POINT(1) TO POINT(NPOINT), WHERE
C     NADRES = (PX-7)*(PY-7) IS THE NUMBER OF PANELS, AND
C     NPOINT .GE. NADRES + M
C
C     .. SCALAR ARGUMENTS ..
      INTEGER IFAIL, M, NADRES, NPOINT, PX, PY
C     .. ARRAY ARGUMENTS ..
      REAL*4 LAMDA(PX), MU(PY), X(M), Y(M)
      INTEGER ADRES(NADRES), POINT(NPOINT)
C     ..
C     .. LOCAL SCALARS ..
C$P 1
      DOUBLE PRECISION SRNAME
      INTEGER ITEMP, K, PANEL, R
C     .. FUNCTION REFERENCES ..
      INTEGER E02ZAZ, P01AAF
C     ..
      DATA SRNAME /8H E02ZAF /
      K = 0
      IF (PX.LT.10) GO TO 40
      ITEMP = PX - 4
      DO 20 R=6,ITEMP
         IF (LAMDA(R-1).GT.LAMDA(R)) K = 1
   20 CONTINUE
   40 IF (PY.LT.10) GO TO 80
      ITEMP = PY - 4
      DO 60 R=6,ITEMP
         IF (MU(R-1).GT.MU(R)) K = 1
   60 CONTINUE
      IF (K.GT.0) GO TO 140
   80 R = (PX-7)*(PY-7)
      IF (M.LE.0 .OR. PX.LT.8 .OR. PY.LT.8 .OR. NADRES.NE.R .OR.NPOINT.
     *LT.R+M) K = 2
      IF (K.GT.0) GO TO 140
      DO 100 PANEL=1,NADRES
         ITEMP = PANEL + M
         ADRES(PANEL) = ITEMP
         POINT(ITEMP) = 0
  100 CONTINUE
      DO 120 R=1,M
         PANEL = (E02ZAZ(PX,LAMDA,X(R))-1)*(PY-7) +E02ZAZ(PY,MU,Y(R))
         ITEMP = ADRES(PANEL)
         ADRES(PANEL) = R
         POINT(ITEMP) = R
         POINT(R) = 0
  120 CONTINUE
      IFAIL = K
      RETURN
  140 IFAIL = P01AAF(IFAIL,K,SRNAME)
      RETURN
      END
      INTEGER FUNCTION E02ZAZ(N, X, T)
C     MARK 6 RELEASE. NAG COPYRIGHT 1977
C
C     GIVEN REAL T AND A REAL ARRAY X OF DIMENSION
C     AT LEAST N - 4, X(5) TO X(N-4) BEING IN NON-
C     DECREASING ORDER, THIS FUNCTION DETERMINES THE
C     INTERVAL X(I+3) TO X(I+4) WHICH CONTAINS T.
C     SPECIFICALLY, E02ZAZ IS ASSIGNED A VALUE I SUCH
C     THAT X(I+3) .LE. T .LT. X(I+4), UNLESS T .LT. X(5)
C     WHEN E02ZAZ = 1, OR T .GE. X(N-4) WHEN E02ZAZ = N-7.
C
C     .. SCALAR ARGUMENTS ..
      REAL*4 T
      INTEGER N
C     .. ARRAY ARGUMENTS ..
      REAL*4 X(N)
C     ..
C     .. LOCAL SCALARS ..
      INTEGER I, L, U
C     ..
      L = 4
      U = N - 3
   20 I = (L+U)/2
      IF (U-L.LE.1) GO TO 60
      IF (T.GE.X(I)) GO TO 40
      U = I
      GO TO 20
   40 L = I
      GO TO 20
   60 E02ZAZ = U - 4
      RETURN
      END

      SUBROUTINE EA06C(A,VALUE,VECTOR,M,IA,IV,W)
C  STANDARD FORTRAN 66 (A VERIFIED PFORT SUBROUTINE)
      REAL A(IA,M),VALUE(M),VECTOR(IV,M),W(1)
      M1=M+1
      M2=M1+M
      W(1)=A(1,1)
      IF(M-2)60,10,15
10    W(2)=A(2,2)
      W(4)=A(2,1)
      GO TO 60
15    CALL MC04B(A,W,W(M1),M,IA,W(M2))
60    CALL EA08C(W,W(M1),VALUE,VECTOR,M,IV,W(M2))
      IF(M.LE.2)RETURN
      DO 56 L=1,M
      DO 56 II=3,M
      I=M-II+1
      M3=M1+I
      IF(W(M3))57,56,57
57    PP=0.
      I1=I+1
      DO 58 K=I1,M
58    PP=PP+A(I,K)*VECTOR(K,L)
      PP=PP/(A(I,I+1)*W(M3))
      DO 59 K=I1,M
59    VECTOR(K,L)=VECTOR(K,L)+PP*A(I,K)
56    CONTINUE
      RETURN
      END
C
      SUBROUTINE EA08C(A,B,VALUE,VEC,M,IV,W)
C  STANDARD FORTRAN 66 (A VERIFIED PFORT SUBROUTINE)
      DIMENSION A(M),B(M),VALUE(M),VEC(1),W(1)
      DATA EPS/1.E-6/,A34/0.0/
C     THIS USES QR ITERATION TO FIND THE EIGENVALUES AND EIGENVECTORS
C  OF THE SYMMETRIC TRIDIAGONAL MATRIX WHOSE DIAGONAL ELEMENTS ARE
C  A(I),I=1,M AND OFF-DIAGONAL ELEMENTS ARE B(I),I=2,M.  THE ARRAY
C  W IS USED FOR WORKSPACE AND MUST HAVE DIMENSION AT LEAST 2*M.
C  WE TREAT VEC AS IF IT HAD DIMENSIONS (IV,M).
      SML=EPS*FLOAT(M)
      CALL EA09C(A,B,W(M+1),M,W)
C     SET VEC TO THE IDENTITY MATRIX.
      DO 5 I=1,M
      VALUE(I)=A(I)
      W(I)=B(I)
      K=(I-1)*IV+1
      L=K+M-1
      DO 3 J=K,L
3     VEC(J)=0.
      KI=K+I-1
5     VEC(KI)=1.
      K=0
      IF(M.EQ.1)RETURN
      DO 200 N3=2,M
      N2=M+2-N3
C     EACH QR ITERATION IS PERFORMED OF ROWS AND COLUMNS N1 TO N2
      MN2=M+N2
      ROOT=W(MN2)
      DO 190 ITER=1,20
      BB=(VALUE(N2)-VALUE(N2-1))*0.5
      CC=W(N2)*W(N2)
      A22=VALUE(N2)
      IF(CC.NE.0.0)A22=A22+CC/(BB+SIGN(1.0,BB)*SQRT(BB*BB+CC))
      DO 125 I=1,N2
      MI=M+I
      IF(ABS(ROOT-A22).LE.ABS(W(MI)-A22))GO TO 125
      ROOT=W(MI)
      MN=M+N2
      W(MI)=W(MN)
      W(MN)=ROOT
125   CONTINUE
      DO 130 II=2,N2
      N1=2+N2-II
      IF(ABS(W(N1)).LE.(ABS(VALUE(N1-1))+ABS(VALUE(N1)))*SML)GO TO 140
130   CONTINUE
      N1=1
140   IF(N2.EQ.N1) GO TO 200
      N2M1=N2-1
      IF(ITER.GE.3)ROOT=A22
      K=K+1
      A22=VALUE(N1)
      A12=A22-ROOT
      A23=W(N1+1)
      A13=A23
      DO180 I=N1,N2M1
      A33=VALUE(I+1)
      IF(I.NE.N2M1)A34=W  (I+2)
      S=SIGN(SQRT(A12*A12+A13*A13),A12)
      SI=A13/S
      CO=A12/S
      JK=I*IV+1
      J1=JK-IV
      J2=J1+MIN0(M,I+K)-1
      DO 160 JI=J1,J2
      V1=VEC(JI)
      V2=VEC(JK)
      VEC(JI)=V1*CO+V2*SI
      VEC(JK)=V2*CO-V1*SI
160   JK=JK+1
      IF(I.NE.N1)  W(I)=S
      A11=CO*A22+SI*A23
      A12=CO*A23+SI*A33
      A13=SI*A34
      A21=CO*A23-SI*A22
      A22=CO*A33-SI*A23
      A23=CO*A34
      VALUE(I)=A11*CO+A12*SI
      A12=-A11*SI+A12*CO
      W(I+1)=A12
180   A22=A22*CO-A21*SI
  190 VALUE(N2)=A22
      WRITE(6,195)
  195 FORMAT(48H1CYCLE DETECTED IN SUBROUTINE EA08 -STOPPING NOW)
       STOP
  200 CONTINUE
C     RAYLEIGH QUOTIENT
      DO 220 J=1,M
      K=(J-1)*IV
      XX=VEC(K+1)**2
      XAX=XX*A(1)
      DO 210 I=2,M
      KI=K+I
      XX=XX+VEC(KI)**2
210   XAX=XAX+VEC(KI)*(2.*B(I)*VEC(KI-1)+A(I)*VEC(KI))
220   VALUE(J)=XAX/XX
      RETURN
      END
C
      SUBROUTINE EA09C(A,B,VALUE,M,OFF)
C  STANDARD FORTRAN 66 (A VERFIED PFORT SUBROUTINE)
      DIMENSION A(M),B(M),VALUE(M),OFF(M)
      DATA A34/0.0/,EPS/1.0E-6/
      SML=EPS*FLOAT(M)
      VALUE(1)=A(1)
      IF(M.EQ.1)RETURN
      DO 10 I=2,M
      VALUE(I)=A(I)
10    OFF(I)=B(I)
C     EACH QR ITERATION IS PERFORMED OF ROWS AND COLUMNS N1 TO N2
      MAXIT=10*M
      DO 90 ITER=1,MAXIT
      DO 45 N3=2,M
      N2=M+2-N3
      DO 30 II=2,N2
      N1=2+N2-II
      IF(ABS(OFF(N1)).LE.(ABS(VALUE(N1-1))+ABS(VALUE(N1)))*SML)GO TO 40
30    CONTINUE
      N1=1
40    IF(N2.NE.N1) GO TO 50
   45 CONTINUE
      RETURN
C     ROOT  IS THE EIGENVALUE OF THE BOTTOM 2*2 MATRIX THAT IS NEAREST
C     TO THE LAST MATRIX ELEMENT AND IS USED TO ACCELERATE THE
C     CONVERGENCE
50    BB=(VALUE(N2)-VALUE(N2-1))*0.5
      CC=OFF(N2)*OFF(N2)
      SBB=1.
      IF(BB.LT.0.)SBB=-1.
      ROOT=VALUE(N2)+CC/(BB+SBB*SQRT(BB*BB+CC))
      N2M1=N2-1
75    A22=VALUE(N1)
      A12=A22-ROOT
      A23=OFF(N1+1)
      A13=A23
      DO 80 I=N1,N2M1
      A33=VALUE(I+1)
      IF(I.NE.N2M1)A34=OFF(I+2)
      S=SQRT(A12*A12+A13*A13)
      SI=A13/S
      CO=A12/S
      IF(I.NE.N1)OFF(I)=S
      A11=CO*A22+SI*A23
      A12=CO*A23+SI*A33
      A13=SI*A34
      A21=CO*A23-SI*A22
      A22=CO*A33-SI*A23
      A23=CO*A34
      VALUE(I)=A11*CO+A12*SI
      A12=-A11*SI+A12*CO
      OFF(I+1)=A12
80    A22=A22*CO-A21*SI
   90 VALUE(N2)=A22
       WRITE(6,100)
  100 FORMAT(39H1LOOPING DETECTED IN EA09-STOPPING NOW )
      STOP
      END
C  VAX LOOKALIKE FOR SUBROUTINE FDATE ON ALLIANT, IRIS, OR GENERAL UNIX.
C  16.1.92 RAC
C
c      	SUBROUTINE FDATE(DAT)
c	CHARACTER  DAT*24,DATVAX*9,TIMVAX*8
c      	DAT='                        '
c	CALL DATE(DATVAX)
c	CALL TIME(TIMVAX)
c      	DAT(7:15) = DATVAX(1:9)
c      	DAT(17:24)= TIMVAX(1:8)
c     	RETURN
c      	END
C
C  Changes to make it work on Alliant 
C
C  CHANGES PUT IN TO MAKE IT WORK ON VAX (THIS VERSION HAS SCALING)
C  1.  12 STATEMENTS FOR CALCULATION OF OVER/UNDERFLOW OF DETERMINANT IN MA21
C      REPLACED BY A SIMPLE ONE WHICH WILL OVERFLOW MORE EASILY(ENTRY MA21CD)
C  2.  CHANGES TO MC10AD TO REPLACE 370-SPECIFIC PARTS....
C      A. U=FLOATI  (6 TIMES)
C      B. NEW ALOG16 PROCEDURE (TWICE)
C      C. SIMPLER 16**DIAG STATEMENT (ONCE)
C  3.  ALL DOUBLE PRECISION REPLACED BY REAL*8 (NOT REALLY NECESSARY).
C  4.  REPLACE A(N),B(N) BY A(1),B(1) IN FM02AD TO AVOID VAX ARRAY CHECKING.
C
      FUNCTION FA01AS(I)
      COMMON/FA01ES/G
      REAL*8 G
      G=DMOD(G* 9228907.D0,4294967296.D0)
      IF(I.GE.0)FA01AS=G/4294967296.D0
      IF(I.LT.0)FA01AS=2.D0*G/4294967296.D0-1.D0
      RETURN
      END
      SUBROUTINE FA01BS(MAX,NRAND)
      NRAND=INT(FA01AS(1)*FLOAT(MAX))+1
      RETURN
      END
      SUBROUTINE FA01CS(IL,IR)
      COMMON/FA01ES/G
      REAL*8 G
      IL=G/65536.D0
      IR=G-65536.D0*FLOAT(IL)
      RETURN
      END
      SUBROUTINE FA01DS(IL,IR)
      COMMON/FA01ES/G
      REAL*8 G
      G=65536.D0*FLOAT(IL)+FLOAT(IR)
      RETURN
      END
C combined these two BLOCK DATA statements so that this file has only one.
      BLOCK DATA
      COMMON/FA01ES/G
      COMMON /MA21ED/LP,JSCALE,EA,EB
C JSCALE NORMALLY =1, TEMPORARILY SET TO 0
      INTEGER * 4 LP/6/,JSCALE/1/
      REAL * 8 EA/1.0D-16/,EB/1.0D-16/
      REAL*8 G
      DATA G/1431655765.D0/
      END
C    FM02AD - A ROUTINE TO COMPUTE THE INNER PRODUCT OF TWO
C      DOUBLE PRECISION REAL VECTORS ACCUMULATING THE RESULT
C      DOUBLE PRECISION.  IT CAN BE USED AS AN ALTERNATIVE
C      TO THE ASSEMBLER VERSION, BUT NOTE THAT IT IS LIKELY
C      TO BE SIGNIFICANTLY SLOWER IN EXECUTION.
C
      REAL*8 FUNCTION FM02AD(N,A,IA,B,IB)
      REAL*8 R1,A,B
C  THE FOLLOWING STATEMENT CHANGED FROM A(N),B(N) TO AVOID VAX DYNAMIC
C  ARRAY CHECK FAILURE.
      DIMENSION A(1),B(1)
C
C    N   THE LENGTH OF THE VECTORS (IF N<= 0  FM02AD = 0)
C    A   THE FIRST VECTOR
C    IA  SUBSCRIPT DISPLACEMENT BETWEEN ELEMENTS OF A
C    B   THE SECOND VECTOR
C    IB  SUBSCRIPT DISPLACEMENT BETWEEN ELEMENTS OF B
C    FM02AD  THE RESULT
C
      R1=0D0
      IF(N.LE.0) GO TO 2
      JA=1
      IF(IA.LT.0) JA=1-(N-1)*IA
      JB=1
      IF(IB.LT.0) JB=1-(N-1)*IB
      I=0
    1 I=I+1
      R1=R1+A(JA)*B(JB)
      JA=JA+IA
      JB=JB+IB
      IF(I.LT.N) GO TO 1
    2 FM02AD=R1
      RETURN
      END
C      BLOCK DATA
C      COMMON /MA21ED/LP,JSCALE,EA,EB
C JSCALE NORMALLY =1, TEMPORARILY SET TO 0
C      INTEGER * 4 LP/6/,JSCALE/1/
C      REAL * 8 EA/1.0D-16/,EB/1.0D-16/
C      END
      SUBROUTINE MA21AD (A,IA,N,B,W,E)
      COMMON /MA21ED/ LP,JSCALE,EA,EB
C*** next statement changed as A and B go to N+1 jms 30.09.98
C***  REAL *  8 A(IA,N),B(N),W(N,N),AA,AC,DET,WW,Q(2),PCK
      REAL *  8 A(IA,N+1),B(N+1),W(N,N),AA,AC,DET,WW,Q(2),PCK
      REAL *8 EPS4/4.0D-16/,ZERO/0.0D0/,ONE/1.0D0/,TWO/2.0D0/,P5/0.5D0/
     1,E,EA,EB,EPSN,XNORM,AXNORM,ENORM,ENORMA,ANORMA,ANORM,ERR,AB,AM
     2,EPS /1.0D-16/,R(4),RA,FM02AD
CTSH++ for G77
      REAL*8 DICT(8)
      CHARACTER*8 TMPDICT(8)/
     1 ' MA21AD ',' MA21BD ',' MA21CD ',' MA21DD ',
     1 ' N IS   ',' PIVOT  ','IS SMALL','IS ZERO '/
      EQUIVALENCE (TMPDICT,DICT)
CTSH--
C     EPS=MACHINE PRECISION,EPS4=EPS*4
C  I3, I4 previously initialised to 0 but changed because of EQUIVALENCE stmt.
C  MRO'D 23/10/81
      INTEGER*4 I3,I4,I1(4)
      REAL*4 UPCK(4)
      LOGICAL*1 L1,L2,L3(4),L4(4),LQ,LA
C      DATA LQ/Z40/	! not needed in this implementation
      EQUIVALENCE (Q(1),R(1),L1),(RA,LA),(WW,II),(UPCK(1),PCK),
     1            (Q(2),L2),(L3(1),I3),(L4(1),I4)
C  Next two lines inserted to initialise I3,I4 by assignment - see above - MRO'D
      I3 = 0
      I4 = 0
      IENT = 1
      GO TO 100
      ENTRY MA21BD(A,IA,N,W,E)
      IENT = 2
      GO TO 100
      ENTRY MA21CD(A,IA,N,DET,IDET,W)
      IENT = 3
      DET = ONE
      IDET = 0
  100 IF(N.LT.1)GO TO 810
      IP = 0
      IF(JSCALE .LE. 0)GO TO 120
C     FIND SCALING FACTORS.
      CALL MC10AD(A,N,IA,W(1,5),W(1,1),IT)
      DO 110 I=1,N
      PCK = W(I,1)
      W(2*I-1,3) = UPCK(1)
  110 W(2*I,3) = UPCK(2)
C     STORE IN W(J,1),J=1,N THE MAXIMAL ELEMENTS  OF COLUMNS AFTER
C     APPLICATION OF SCALING.
  120 EPSN = ZERO
      DO 140 J=1,N
      AB = ZERO
      DO 130 I=1,N
      AM =  DABS(A(I,J))
      IF(JSCALE.GT.0)AM=AM*W(I,3)
  130 AB = DMAX1(AB,AM)
      W(J,1) = AB
      IF(JSCALE.GT.0)AB=AB*W(J,4)
  140 EPSN = DMAX1(EPSN,AB)
      EPSN = EPSN*EPS4
      IF((E.LE.ZERO).OR.(IENT.EQ.3))GO TO 160
C     MAKE COPY OF MATRIX
      DO 150 I=1,N
      DO 150 J=1,N
  150 W(I,J+5) = A(I,J)
C
C     FACTORISATION OF MATRIX INTO L*U, WHERE L IS A LOWER UNIT
C     TRIANGLE AND U IS UPPER TRIANGLE
  160 DO 230 L=1,N
      AM = ZERO
      II = L
C     EVALUATE  ELEMENTS IN PIVOTAL COLUMN.
      DO 170 I=L,N
      A(I,L)=A(I,L)-FM02AD(L-1,A(I,1),IA,A(1,L),1)
C     LOOK FOR MAXIMUM ELEMENT IN PIVOTAL COLUMN.
      AB =  DABS(A(I,L))
      IF(JSCALE .GT. 0)AB = AB*W(I,3)
      IF(AM .GE. AB)GO TO 170
      AM = AB
      II = I
  170 CONTINUE
C
C     TEST FOR SMALL OR ZERO PIVOT.
      AB = W(L,1)
      IF(AM .LE. EPS4*AB)IP=-L
      IF(AM.NE.ZERO)GO TO 180
      IF(IENT.EQ.2)GO TO 820
      IP=L
  180 IF(IENT.NE.3)GO TO 190
C
C     THE NEXT 12 STATEMENTS CALCULATE THE DETERMINANT OF MATRIX A.
C     TO PREVENT OVERFLOWS/UNDERFLOWS THE EXPONENTS OF THE NUMBERS
C     BEING MULTIPLIED ARE EXAMINED AND THE EXCESS EXPONENT IS STORED
C     IN IDET.
C FOR VAX NEXT TWELVE STATEMENTS REPLACED- MAY NEED ATTENTION IF MA21CD
C IS EVER NEEDED FOR BIG MATRICES.
C      Q(1)= DABS(DET)
C      Q(2)= DABS(A(II,L))
C      L3(4) = L1
C      L4(4) = L2
C      K = IDET+I3+I4-128
C      I3 = K
C      L2 = LQ
C      IF(IABS(K) .GT. 62)K=ISIGN(62,I3)
C      IDET = I3-K
C      I3 = K+64
C      L1 = L3(4)
C      DET =DSIGN(Q(1),DET)*DSIGN(Q(2),A(II,L))
      DET=DET*A(II,L)
C
  190 IF(II .EQ. L)GO TO 220
C
C     INTERCHANGE ROWS N AND II
C     INTERCHANGE EQUILIBRATION FACTORS
C     W(L,1)=II MEANS INTERCHANGE BETWEEN ROWS L AND II
      IF(IENT.EQ.3)DET=-DET
      IF(JSCALE .LE. 0)GO TO 200
      AA = W(L,3)
      W(L,3) = W(II,3)
      W(II,3) = AA
  200 DO 210 I=1,N
      AA = A(L,I)
      A(L,I) = A(II,I)
  210 A(II,I) = AA
  220 W(L,1) = WW
      IF( L .EQ. N)GO TO 240
      AA = ONE
      AC = A(L,L)
      IF( DABS(AC) .NE. ZERO)AA = AA/AC
      K= L+1
C     UPPER TRIANGLE
      DO 230 I=K,N
      A(I,L) = A(I,L)*AA
  230 A(L,I)=A(L,I)-FM02AD(L-1,A(L,1),IA,A(1,I),1)
C
  240 IF(IENT -2)250,500,720
C
  250 IF(E .LE. ZERO)GO TO 270
      DO 260 I=1,N
      W(I,5) = B(I)
  260 W(I,2) = ZERO
      IT = 0
C
C     FORWARD SUBSTITUTION
  270 DO 280 I=1,N
      WW = W(I,1)
      AA = B(II)
      B(II) = B(I)
  280 B(I)=AA-FM02AD(I-1,A(I,1),IA,B(1),1)
      ENORMA = ENORM
      ENORM = ZERO
C
C     CALCULATE NORMS OF X,CHANGE IN SOLUTION
C     BACKWARD SUBSTITUTION
      DO 310 K=1,N
      I=N+1-K
      AA = A(I,I)
      IF(   DABS(AA) .EQ. ZERO)GO TO 290
C*** This statement gives problems as i+1 > 2nd dimension
      AC=B(I)-FM02AD(N-I,A(I,I+1),IA,B(I+1),1)
      B(I) = AC/AA
      GO TO 300
  290 IP=I
      B(I) = ZERO
  300 AM=ONE
      IF(JSCALE.GT.0)AM=W(I,4)
  310 ENORM=DMAX1(ENORM, DABS(B(I))/AM)
      IF((E .LE. ZERO) .OR. (IP .NE. 0 ))GO TO 720
      IF(IT .EQ. 0)GO TO 320
      IF(ENORM .GT. P5*ENORMA) GO TO 460
      IF(ENORM   -  EPS*XNORM)470,470,330
  320 XNORM=ENORM
      CALL FA01CS(IRANDL,IRANDR)
      CALL FA01DS(21845,21845)
C
C     UPDATE SOLUTION VECTOR X
  330 DO 340 I=1,N
  340 W(I,2) = W(I,2)+B(I)
      IT = IT+1
C
C     COMPUTE RESIDUAL
  350 DO 450 J=1,N
      IF(IENT .NE. 2)GO TO 360
      AA = ZERO
      IF(L .EQ. J)AA = ONE
      GO TO 370
  360 AA = W(J,5)
  370 AC=AA-FM02AD(N,W(J,6),N,W(1,2),1)
      AA = ZERO
      IF(EA)400,430,380
C
C     MAKE PSEUDO RANDOM CHANGES TO ELEMENTS OF A AND B
  380 DO 390 K=1,N
  390 AA = AA+FA01AS(-K)*W(J,K+5)*W(K,2)
      GO TO 420
  400 DO 410 K=1,N
  410 AA = AA+FA01AS(-K)*W(K,2)
  420 AC =AC-DABS(EA)*AA
  430 IF(IENT .EQ. 2)GO TO 440
      AA  = DABS(EB)*FA01AS(-J)
      IF(EB .GE. ZERO)AA=AA*W(J,5)
      AC = AA+AC
      B(J) = AC
      GO TO 450
  440 W(J,5) = AC
  450 CONTINUE
      IF(IENT-2)270,630,270
  460 ENORM = ENORMA
  470 DO 480 I=1,N
  480 B(I) = W(I,2)
      E=ENORM
      IF(JSCALE .LE. 0)GO TO 710
C
C     SET UP ACCURACY ESTIMATES FOR SOLUTION VECTOR.
      ERR = ZERO
      DO 490 I=1,N
      W(I,2) = ENORM*W(I,4)
      AB = W(I,2)
  490 ERR = DMAX1(ERR,AB)
      E=ERR
      GO TO 710
C
C     OVERWRITE LU FACTORISATION OF A BY INVERSE OF PERMUTED A.
  500 IF(N .LT. 2)GO TO 520
      DO 510 I=2,N
      K=I-1
      DO 510 J=1,K
  510 A(I,J)=-A(I,J)-FM02AD(I-1-J,A(I,J+1),IA,A(J+1,J),1)
  520 DO 540 K=1,N
      I = N+1-K
      ERR = ONE/EPSN
      IF(JSCALE .GT. 0)ERR = ERR*W(I,4)
      DO 530 J=I,N
      W(J,2) = A(I,J)
  530 A(I,J) = ZERO
      A(I,I) = ONE
      AA = ONE/W(I,2)
      DO 540 J=1,N
      AB = ONE
      IF (JSCALE .GT. 0)AB = W(J,3)
      AC=A(I,J)-FM02AD(N-I,W(I+1,2),1,A(I+1,J),1)
      A(I,J) = AC*AA
      IF( DABS(A(I,J)) .GE. ERR*AB)IP=N*I+J
  540 CONTINUE
      ANORM = ZERO
      DO 590 K=1,N
      I = N+1-K
      WW = W(I,1)
      IF(II .EQ. I)GO TO 570
      IF(JSCALE .LE. 0)GO TO 550
      AA = W(I,3)
      W(I,3) = W(II,3)
      W(II,3) = AA
  550 DO 560 J=1,N
      AA = A(J,II)
      A(J,II) = A(J,I)
  560 A(J,I) = AA
      W(I,1) = W(II,1)
  570 ENORM = ZERO
      DO 580 J=1,N
      AB =  DABS(A(J,II))
      IF(JSCALE.GT.0)AB=AB/W(J,4)
  580 ENORM = DMAX1(ENORM,AB)
      W(II,1) = ENORM
      AB = ONE
      IF(JSCALE.GT.0)AB=W(II,3)
  590 ANORM = DMAX1(ANORM,ENORM/AB)
      IF((E .LE. ZERO) .OR. (IP .NE. 0 ))GO TO 720
       CALL FA01CS(IRANDL,IRANDR)
      CALL FA01DS(21845,21845)
      AXNORM=ANORM
  600 ANORMA = ANORM
      ANORM = ZERO
      L=0
  610 L=L+1
C
C     INVERSE OF A IS ITERATIVELY REFINED BY COLUMNS
C     MAKE COPY OF APPROPIATE COLUMN.
      DO 620 I=1,N
  620 W(I,2) = A(I,L)
      GO TO 350
C
C     CALCULATE THE CHANGE IN APPROPIATE COLUMN OF INVERSE OF A.
  630 ENORM = ZERO
      DO 640 I=1,N
      W(I,2) = ZERO
      W(I,2)=W(I,2)+FM02AD(N,A(I,1),IA,W(1,5),1)
      AM=ONE
      IF(JSCALE.GT.0)AM=W(I,4)
  640 ENORM=DMAX1(ENORM, DABS(W(I,2))/AM)
      AB = W(L,1)
      IF(ENORM .GT. P5*AB)GO TO 660
C
C     UPDATE APPROPIATE COLUMN OF INVERSE OF A.
      DO 650 J=1,N
  650 A(J,L) = A(J,L)+W(J,2)
      W(L,1) = ENORM
      AB = ENORM
  660 ERR = ONE
      IF(JSCALE .GT. 0)ERR = W(L,3)
      ANORM = DMAX1(ANORM,AB/ERR)
      IF(L.LT.N)GO TO 610
      IF(ANORM .GT. P5*ANORMA)GO TO 670
      IF(ANORM   -  EPS*AXNORM)680,680,600
  670 ANORM = ANORMA
  680 IF(JSCALE .LE. 0)GO TO 700
C
C     SET UP ACCURACY ESTIMATES FOR MATRIX INVERSE.
      ENORM = ZERO
      ERR = ZERO
      DO 690 I=1,N
      AB = ANORM*W(I,3)
      ENORM = DMAX1(ENORM,AB)
      W(I,1) = AB
      AB = W(I,4)
      ERR = DMAX1(ERR,AB)
  690 W(I,2) = AB
      E=ENORM*ERR
      GO TO 710
  700 E=ANORM
  710 CALL FA01DS(IRANDL,IRANDR)
      GO TO 850
C
      ENTRY MA21DD (A,IA,N,B,W,E)
      IENT = 4
      IF ( N .LE. 0)GO TO 810
C
C     CHECK WHETHER ON A PREVIOUS ENTRY A SMALL PIVOT WAS FOUND,IF
C     SO PUT ON ERROR FLAG.
      IP = 0
      WW = W(N,1)
      IF(II .GT. 0)GO TO 250
      IP = II
      I = -II
      II = N
      W(N,1) = WW
      GO TO 250
C
C     THE REMAINING INSTRUCTIONS HANDLE ERROR DIAGNOSTICS,ETC.
  720 IF(IP)730,800,740
  730 II=IP
      W(N,1) = WW
      J=7
      GO TO 770
  740 IF(IP.LE.N)GO TO 760
      I=(IP-1)/N
      J=IP-N*I
      WRITE(LP,750)DICT(IENT),I,J
  750 FORMAT('::',A8,'HAS FOUND THAT INVERSE ELEMENT (',I4,',',I4,') IS LARGE
     1,RESULTS MAY BE UNRELIABLE')
      GO TO 790
  760 J=8
  770 I=IABS(IP)
      WRITE(LP,780)DICT(IENT),DICT(6),I,DICT(J)
  780 FORMAT('::',A8,'HAS FOUND THAT',A8,I4,1X,A8,',RESULTS MAY',
     1  ' BE UNRELIABLE')
  790 E=-ONE
      GO TO 850
  800 E=ZERO
      GO TO 850
  810 WRITE(LP,830)DICT(IENT),DICT(5),N
      GO TO 840
  820 WRITE(LP,830)DICT(IENT),DICT(6),L,DICT(8 )
  830 FORMAT(':: ERROR RETURN FROM',2A8,I5,1X,A8)
  840 E=-TWO
  850 RETURN
      END
      SUBROUTINE MC10AD(A,N,NN,DIAG,RES,IS)
      REAL*8 A(NN,NN)
      REAL*4 RES(N,4)
      INTEGER*2DIAG(N,2),JU(2),KU
C     RES IS USED TO RETURN SCALING FACTORS AS INTEGRAL
C          POWERS OF BASE, AND AS WORKSPACE
C     IS IS SET TO 0 ON SUCCESSFUL COMPLETION, TO I IF ROW I HAS ONLY
C        ZERO ELEMENTS, TO -I IF COLUMN I HAS ONLY ZERO ELEMENTS
C     DIAG IS USED TO HOLD COUNTS OF NON-ZEROS IN ROWS AND COLUMNS
C      AND TO RETURN SCALING POWERS
C
      DATA SMIN/.01/
C     SMIN IS USED IN A CONVERGENCE TEST ON (RESIDUAL NORM)**2
      LOGICAL*1 IU,IW(3)
C      EQUIVALENCE (UU,IW(1),KU),(U,IU,JU(1)) ......370 SPECIFIC
C     SET UP CONSTANTS
C      UU = 100.  ......370 SPECIFIC
      UU=ALOG(16.0)
      IS=0
C     INITIALISE FOR ACCUMULATION OF SUMS AND PRODUCTS
      DO 2 L=1,2
      DO 2 I=1,N
      RES(I,L)=0.
      RES(I,L+2)=0.
    2 DIAG(I,L)=0
      DO 3 I=1,N
      DO 3 J=1,N
      U=DABS(A(I,J))
      IF(U.EQ.0.)GO TO 3
C     ON THE IBM 360 THE FOLLOWING TWO INSTRUCTIONS FIND THE SMALLEST
C     INTEGER GREATER THAN ALOG16(U).
C      IW(2) = IU
C      U=UU-64.
      U=ALOG(U)/UU
      U=AINT(U+1.0)
C     COUNT NON-ZEROS IN ROW AND COLUMN
      DIAG(I,1)=DIAG(I,1)+1
      DIAG(J,2)=DIAG(J,2)+1
      RES(I,1)=RES(I,1)+U
      RES(J,3)=RES(J,3)+U
    3 CONTINUE
C     COMPUTE RHS VECTORS TESTING FOR ZERO ROW OR COLUMN
      SSUM=0.
      J=0
C      JU(1)=17920  .......370 SPECIFIC
      DO 8 I=1,N
      J=J+DIAG(I,1)
      DO 9 L=1,2
      IF(DIAG(I,L).GT.0 )GO TO 153
      DIAG(I,L)=1
      IS=I*(3-2*L)
  153 CONTINUE
C      ON IBM 360 NEXT INSTRUCTION SETS U TO VALUE OF POSITIVE INTEGER
C      JU(2)=DIAG(I,L)
CTSH      U=FLOATI(DIAG(I,L))
CTSH++
      U=DIAG(I,L)
CTSH--
    9 RES(I,2*L-1)=RES(I,2*L-1)/U
    8 SSUM=SSUM+RES(I,3)
      SM=SMIN*J
C     SWEEP TO COMPUTE INITIAL RESIDUAL VECTOR
      RSUM=0.
      DO 110 I=1,N
      SUM = SSUM
      IF(DIAG(I,1).GE. N)GO TO 109
      SUM=0.
      DO 10 J=1,N
      IF(A(I,J).EQ.0D0)GO TO 10
      SUM=SUM+RES(J,3)
   10 CONTINUE
C      ON IBM 360 NEXT INSTRUCTION SETS U TO VALUE OF POSITIVE INTEGER
C  109 JU(2)=DIAG(I,1)
CTSH  109 U=FLOATI(DIAG(I,1))
CTSH++
  109 U=DIAG(I,1)
CTSH--
      RES(I,1)=RES(I,1)-SUM/U
  110 RSUM=RSUM+RES(I,1)
C     INITIALISE ITERATION
      E=0.
      E1=0.
      Q=1.
      S=0.
      DO 11 I=1,N
C      ON IBM 360 NEXT INSTRUCTION SETS U TO VALUE OF POSITIVE INTEGER
C      JU(2)=DIAG(I,1)
CTSH      U=FLOATI(DIAG(I,1))
CTSH++
      U=DIAG(I,1)
CTSH--
   11 S=S+U*RES(I,1)**2
      L=2
      IF(S.LE.SM)GO TO 100
C     ITERATION STEP
   20 EM=E*E1
C    SWEEP THROUGH MATRIX TO UPDATE RESIDUAL VECTOR
      DO 220 I=1,N
      SUM=RSUM
      IF(DIAG(I,L).GE. N)GO TO 220
      SUM=0.
      DO 22 J=1,N
      IF(L.EQ.2)GO TO 21
      IF(A(I,J))19,22,19
   21 IF(A(J,I))19,22,19
   19 SUM=SUM+RES(J,3-L)
   22 CONTINUE
  220 RES(I,L)=RES(I,L)+SUM
      S1=S
      S=0.
      RSUM=0.
      DO 23 I=1,N
      V=-RES(I,L)/Q
C      ON IBM 360 NEXT INSTRUCTION SETS U TO VALUE OF POSITIVE INTEGER
C      JU(2)=DIAG(I,L)
CTSH      U=FLOATI(DIAG(I,L))
CTSH++
      U=DIAG(I,L)
CTSH--
      RES(I,L)=V/U
      RSUM=RSUM+RES(I,L)
   23 S=S+V*RES(I,L)
      E1=E
      E=Q*S/S1
      Q1=Q
      Q=1.-E
      M=3-L
      IF(S.GT.SM)GO TO 27
      E=M-1
      M=1
      Q=1.
   27 IF(L.EQ.2)GO TO 25
      QM=Q*Q1
      DO 24 I=1,N
      RES (I,4)=(EM*RES(I,4)+RES(I,2))/QM
   24 RES(I,3)=RES(I,3)+RES(I,4)
   25 L=M
      DO 26 I=1,N
C      ON IBM 360 NEXT INSTRUCTION SETS U TO VALUE OF POSITIVE INTEGER
C      JU(2)=DIAG(I,L)
CTSH      U=FLOATI(DIAG(I,L))
CTSH++
      U=DIAG(I,L)
CTSH--
   26 RES(I,L)=RES(I,L)*U*E
      IF(S.GT.SM  )GO TO 20
C      SWEEP THROUGH MATRIX TO GET ROW SCALING POWERS
  100 DO 103 I=1,N
      DO 103 J=1,N
      U=DABS(A(I,J))
      IF(U.EQ.0.)GO TO 103
C      ON IBM 360 NEXT TWO INSTRUCTIONS FIND THE SMALLEST INTEGER
C     LESS THAN ALOG16(U)
C      IW(2)=IU
C      U=UU-64.
      U=ALOG(U)/UU
      U=AINT(U+1.0)
      RES(I,1)=RES(I,1)+RES(J,3)-U
  103 CONTINUE
C      CONVERT POWERS TO INTEGERS
C      JU(1)=17920  .......370 SPECIFIC
      DO 104 I=1,N
C      ON IBM 360 NEXT INSTRUCTION SETS U TO VALUE OF POSITIVE INTEGER
C      JU(2)=DIAG(I,1)
CTSH      U=FLOATI(DIAG(I,1))
CTSH++
      U=DIAG(I,1)
CTSH--
      V=RES(I,1)/U
      DIAG(I,1)=V+SIGN(0.5,V)
  104 DIAG(I,2)=-(RES(I,3)+SIGN(0.5,RES(I,3)))
C      SET SCALING FACTORS IN RES
C      U=1.   .......370 SPECIFIC
      DO 105 L=1,2
      DO 105 I=1,N
  105 RES(I,L)=16.0**DIAG(I,L)
C     THE FOLLOWING THREE STATEMENTS WERE EQUIVALENT TO THE ONE REPLACING IT
C      KU=DIAG(I,L)+65
C      IU=IW(2)
C  105 RES(I,L)=U
      RETURN
      END
C
      SUBROUTINE MC04B(A,ALPHA,BETA,M,IA,Q)
C STANDARD FORTRAN 66 (A VERIFIED PFORT SUBROUTINE)
      DIMENSION A(IA,1),ALPHA(1),BETA(1),Q(1)
      ALPHA(1)=A(1,1)
      DO 21 J=2,M
      J1=J-1
      DO 22 I=1,J1
      A(I,J)=A(J,I)
   22 CONTINUE
      ALPHA(J)=A(J,J)
   21 CONTINUE
      M1=M-1
      M2=M-2
      DO 1 I=1,M2
      PP=0.0
      I1=I+1
      DO 2 J=I1,M
      PP=PP+A(I,J)**2
    2 CONTINUE
      PP1=SQRT(PP)
      IF(A(I,I+1))3,5,5
    5 BETA(I+1)=-PP1
      GO TO 6
    3 BETA(I+1)=PP1
    6 IF(PP)1,1,17
   17 H=PP-BETA(I+1)*A(I,I+1)
      A(I,I+1)=A(I,I+1)-BETA(I+1)
      DO 7 KI=I1,M
      QJ=0.0
      DO 8 KJ=I1,KI
      QJ=QJ+A(KJ,KI)*A(I,KJ)
    8 CONTINUE
      IF(KI-M)19,20,20
   19 I2=KI+1
      DO 18 KJ=I2,M
      QJ=QJ+A(KI,KJ)*A(I,KJ)
   18 CONTINUE
   20 Q(KI)=QJ/H
    7 CONTINUE
      BIGK=0.0
      DO 9 KJ=I1,M
      BIGK=BIGK+A(I,KJ)*Q(KJ)
    9 CONTINUE
      BIGK=BIGK/(2.0*H)
      DO 10 KJ=I1,M
      Q(KJ)=Q(KJ)-BIGK*A(I,KJ)
   10 CONTINUE
      DO 11 KI=I1,M
      DO 12 KJ=KI,M
      A(KI,KJ)=A(KI,KJ)-Q(KI)*A(I,KJ)-Q(KJ)*A(I,KI)
   12 CONTINUE
   11 CONTINUE
    1 CONTINUE
      DO 23 I=2,M
      H=ALPHA(I)
      ALPHA(I)=A(I,I)
      A(I,I)=H
   23 CONTINUE
      BETA(M)=A(M-1,M)
      RETURN
      END
      INTEGER FUNCTION P01ABF(IFAIL,IERROR,SRNAME,NREC,REC)
C     MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986.
C     MARK 13 REVISED. IER-621 (APR 1988).
C     MARK 13B REVISED. IER-668 (AUG 1988).
C
C     P01ABF is the error-handling routine for the NAG Library.
C
C     P01ABF either returns the value of IERROR through the routine
C     name (soft failure), or terminates execution of the program
C     (hard failure). Diagnostic messages may be output.
C
C     If IERROR = 0 (successful exit from the calling routine),
C     the value 0 is returned through the routine name, and no
C     message is output
C
C     If IERROR is non-zero (abnormal exit from the calling routine),
C     the action taken depends on the value of IFAIL.
C
C     IFAIL =  1: soft failure, silent exit (i.e. no messages are
C                 output)
C     IFAIL = -1: soft failure, noisy exit (i.e. messages are output)
C     IFAIL =-13: soft failure, noisy exit but standard messages from
C                 P01ABF are suppressed
C     IFAIL =  0: hard failure, noisy exit
C
C     For compatibility with certain routines included before Mark 12
C     P01ABF also allows an alternative specification of IFAIL in which
C     it is regarded as a decimal integer with least significant digits
C     cba. Then
C
C     a = 0: hard failure  a = 1: soft failure
C     b = 0: silent exit   b = 1: noisy exit
C
C     except that hard failure now always implies a noisy exit.
C
C     S.Hammarling, M.P.Hooper and J.J.du Croz, NAG Central Office.
C
C     .. Scalar Arguments ..
      INTEGER                 IERROR, IFAIL, NREC
      CHARACTER*(*)           SRNAME
C     .. Array Arguments ..
      CHARACTER*(*)           REC(*)
C     .. Local Scalars ..
      INTEGER                 I, NERR
      CHARACTER*72            MESS
C     .. External Subroutines ..
      EXTERNAL                P01ABZ, X04AAF, X04BAF
C     .. Intrinsic Functions ..
      INTRINSIC               ABS, MOD
C     .. Executable Statements ..
      IF (IERROR.NE.0) THEN
C        Abnormal exit from calling routine
         IF (IFAIL.EQ.-1 .OR. IFAIL.EQ.0 .OR. IFAIL.EQ.-13 .OR.
     *       (IFAIL.GT.0 .AND. MOD(IFAIL/10,10).NE.0)) THEN
C           Noisy exit
            CALL X04AAF(0,NERR)
            DO 20 I = 1, NREC
               CALL X04BAF(NERR,REC(I))
   20       CONTINUE
            IF (IFAIL.NE.-13) THEN
               WRITE (MESS,FMT=7786) SRNAME, IERROR
               CALL X04BAF(NERR,MESS)
               IF (ABS(MOD(IFAIL,10)).NE.1) THEN
C                 Hard failure
                  CALL X04BAF(NERR,
     *                     ' ** NAG hard failure - execution terminated'
     *                        )
                  CALL P01ABZ
               ELSE
C                 Soft failure
                  CALL X04BAF(NERR,
     *                        ' ** NAG soft failure - control returned')
               END IF
            END IF
         END IF
      END IF
      P01ABF = IERROR
      RETURN
C
 7786 FORMAT (':: ** ABNORMAL EXIT from NAG Library routine ',A,': IFAIL',
     *  ' =',I6)
      END
      SUBROUTINE P01ABZ
C     MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986.
C
C     Terminates execution when a hard failure occurs.
C
C     ******************** IMPLEMENTATION NOTE ********************
C     The following STOP statement may be replaced by a call to an
C     implementation-dependent routine to display a message and/or
C     to abort the program.
C     *************************************************************
C     .. Executable Statements ..
      STOP
      END
* P01FTEXT
*UPTODATE P01AAFTEXT
      INTEGER FUNCTION P01AAF(IFAIL, ERROR, SRNAME)
C     MARK 1 RELEASE.  NAG COPYRIGHT 1971
C     MARK 3 REVISED
C     MARK 4A REVISED, IER-45
C     MARK 4.5 REVISED
C     MARK 7 REVISED (DEC 1978)
C     RETURNS THE VALUE OF ERROR OR TERMINATES THE PROGRAM.
      INTEGER ERROR, IFAIL, NOUT
C$P 1
      DOUBLE PRECISION SRNAME
C     TEST IF NO ERROR DETECTED
      IF (ERROR.EQ.0) GO TO 20
C     DETERMINE OUTPUT UNIT FOR MESSAGE
      CALL X04AAF (0,NOUT)
C     TEST FOR SOFT FAILURE
      IF (MOD(IFAIL,10).EQ.1) GO TO 10
C     HARD FAILURE
      WRITE (NOUT,7786) SRNAME, ERROR
C     STOPPING MECHANISM MAY ALSO DIFFER
C     CALL P01AAZ
      STOP
C     SOFT FAIL
C     TEST IF ERROR MESSAGES SUPPRESSED
   10 IF (MOD(IFAIL/10,10).EQ.0) GO TO 20
      WRITE (NOUT,7786) SRNAME, ERROR
   20 P01AAF = ERROR
      RETURN
 7786 FORMAT (1H0, 38HERROR DETECTED BY NAG LIBRARY ROUTINE , A8,
     * 11H - IFAIL = , I5//)
      END
** END OF P01AAFTEXT
      DOUBLE PRECISION FUNCTION S18AEF(X,IFAIL)
C     MARK 7 RELEASE. NAG COPYRIGHT 1978.
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C     BESSEL FUNCTION I0(X)
C
C     **************************************************************
C
C     TO EXTRACT THE CORRECT CODE FOR A PARTICULAR MACHINE-RANGE,
C     ACTIVATE THE STATEMENTS CONTAINED IN COMMENTS BEGINNING  CDD ,
C     WHERE  DD  IS THE APPROXIMATE NUMBER OF SIGNIFICANT DECIMAL
C     DIGITS REPRESENTED BY THE MACHINE
C     DELETE THE ILLEGAL DUMMY STATEMENTS OF THE FORM
C     * EXPANSION (NNNN) *
C
C     ALSO INSERT APPROPRIATE DATA STATEMENTS TO DEFINE CONSTANTS
C     WHICH DEPEND ON THE RANGE OF NUMBERS REPRESENTED BY THE
C     MACHINE, RATHER THAN THE PRECISION (SUITABLE STATEMENTS FOR
C     SOME MACHINES ARE CONTAINED IN COMMENTS BEGINNING CRD WHERE
C     D IS A DIGIT WHICH SIMPLY DISTINGUISHES A GROUP OF MACHINES).
C     DELETE THE ILLEGAL DUMMY DATA STATEMENTS WITH VALUES WRITTEN
C     *VALUE*
C
C     **************************************************************
C
C     .. Parameters ..
      CHARACTER*6                      SRNAME
      PARAMETER                        (SRNAME='S18AEF')
C     .. Scalar Arguments ..
      DOUBLE PRECISION                 X
      INTEGER                          IFAIL
C     .. Local Scalars ..
      DOUBLE PRECISION                 G, T, XBIG, XVSMAL, Y, YBIG
C     .. Local Arrays ..
      CHARACTER*1                      P01REC(1)
C     .. External Functions ..
      INTEGER                          P01ABF
      EXTERNAL                         P01ABF
C     .. Intrinsic Functions ..
      INTRINSIC                        ABS, LOG, EXP
C     .. Data statements ..
C08   DATA XVSMAL /1.0D-4/
C09   DATA XVSMAL /3.2D-5/
C12   DATA XVSMAL /1.0D-6/
C15   DATA XVSMAL /3.2D-8/
      DATA XVSMAL /3.2D-9/
C19   DATA XVSMAL /3.2D-10/
C     FOR DEC VAX
      DATA XBIG,YBIG / 88.0D0, 1.7D+37 /
C     XBIG = LARGEST X SUCH THAT  EXP(X)/SQRT(X) .LT. MAXREAL
C             (ROUNDED DOWN)
C     YBIG = EXP(XBIG)/SQRT(XBIG)
C     FOR IEEE SINGLE PRECISION
CR0   DATA XBIG,YBIG /90.9E0,3.1E+38/
C     FOR IBM 360/370 AND SIMILAR MACHINES
CR1   DATA XBIG,YBIG /177.2D0,2.7D+75/
C     FOR DEC-10, HONEYWELL, UNIVAC 1100(S.P.)
CR2   DATA XBIG,YBIG /90.2D0,6.3D+37/
C     FOR ICL 1900
CR3   DATA XBIG,YBIG /179.3D0,2.2D+76/
C     FOR CDC 7600/CYBER
CR4   DATA XBIG,YBIG /744.9D0,4.7D+321/
C     FOR UNIVAC 1100(D.P.)
CR5   DATA XBIG,YBIG /712.3D0,3.3D+307/
C     FOR IEEE DOUBLE PRECISION
CR7   DATA XBIG,YBIG /713.0D0,1.7D+308/
C     .. Executable Statements ..
C
      T = ABS(X)
C     ERROR 1 TEST
      IF (T.GT.XBIG) GO TO 80
      IFAIL = 0
C     X RANGE TEST
      IF (T.GT.12.0D0) GO TO 60
C     TEST FOR MIDDLE RANGE
      IF (T.GT.4.0D0) GO TO 40
C     SMALL X
C     TEST FOR VERY SMALL X
      IF (T.GT.XVSMAL) GO TO 20
      S18AEF = 1.0D0
      GO TO 100
   20 G = EXP(T)
      T = 0.5D0*T - 1.0D0
C
C
C     EXPANSION (0035) EVALUATED AS Y(T)  --PRECISION 08E.09
C08   Y = ((((((((((((((+1.25170153D-5)*T-4.67023479D-5)
C08  *    *T+1.19121179D-4)*T-3.77327453D-4)*T+1.16281590D-3)
C08  *    *T-3.14217301D-3)*T+7.69874623D-3)*T-1.71264378D-2)
C08  *    *T+3.41510170D-2)*T-6.04328134D-2)*T+9.41615739D-2)
C08  *    *T-1.28895516D-1)*T+1.57686847D-1)*T-1.86478069D-1)*T +
C08  *    3.08508323D-1
C
C     EXPANSION (0035) EVALUATED AS Y(T)  --PRECISION 09E.10
C09   Y = (((((((((((((((-3.149367911D-6)*T+1.251701532D-5)
C09  *    *T-3.489221823D-5)*T+1.191211786D-4)*T-3.950426472D-4)
C09  *    *T+1.162815900D-3)*T-3.128640567D-3)*T+7.698746228D-3)
C09  *    *T-1.713197375D-2)*T+3.415101696D-2)*T-6.043165089D-2)
C09  *    *T+9.416157393D-2)*T-1.288956234D-1)*T+1.576868468D-1)
C09  *    *T-1.864780666D-1)*T + 3.085083225D-1
C
C     EXPANSION (0035) EVALUATED AS Y(T)  --PRECISION 12E.13
C12   Y = (((((((((((((((+3.552213748272D-8*T-1.672619002648D-7)
C12  *    *T+5.866981988794D-7)*T-2.438504834571D-6)
C12  *    *T+9.830542086049D-6)*T-3.613622861718D-5)
C12  *    *T+1.236706911308D-4)*T-3.938874947045D-4)
C12  *    *T+1.158888468912D-3)*T-3.129251464727D-3)
C12  *    *T+7.700609104316D-3)*T-1.713179048190D-2)
C12  *    *T+3.415053905809D-2)*T-6.043168004160D-2)
C12  *    *T+9.416163400294D-2)*T-1.288956212997D-1)
C12  *    *T+1.576868439705D-1
C12   Y = (Y*T-1.864780666100D-1)*T + 3.085083225537D-1
C
C     EXPANSION (0035) EVALUATED AS Y(T)  --PRECISION 15E.16
C15   Y = (((((((((((((((-2.530899946655355D-10
C15  *    *T+1.378625787804520D-9)*T-5.841016145715136D-9)
C15  *    *T+2.862900854369717D-8)*T-1.361952673926676D-7)
C15  *    *T+6.013460978748421D-7)*T-2.502852637768319D-6)
C15  *    *T+9.813309263701611D-6)*T-3.606463713226221D-5)
C15  *    *T+1.236829425904201D-4)*T-3.939345014889662D-4)
C15  *    *T+1.158883078270167D-3)*T-3.129232874073100D-3)
C15  *    *T+7.700610548238343D-3)*T-1.713179479244525D-2)
C15  *    *T+3.415053883594465D-2)*T-6.043167950086759D-2
C15   Y = ((((Y*T+9.416163402029200D-2)*T-1.288956213305214D-1)
C15  *    *T+1.576868439699908D-1)*T-1.864780666094668D-1)*T +
C15  *    3.085083225536711D-1
C
C     EXPANSION (0035) EVALUATED AS Y(T)  --PRECISION 17E.18
      Y = (((((((((((((((-7.48150165756234957D-12
     *    *T+4.44484446637868974D-11)*T-2.10071360134551962D-10)
     *    *T+1.13415934215369209D-9)*T-5.94856273204259507D-9)
     *    *T+2.92096163521178835D-8)*T-1.36042013507151017D-7)
     *    *T+6.00566861079330132D-7)*T-2.50298975966588680D-6)
     *    *T+9.81395862769787105D-6)*T-3.60645571444886286D-5)
     *    *T+1.23682594989692688D-4)*T-3.93934532072526720D-4)
     *    *T+1.15888319775791686D-3)*T-3.12923286656374358D-3)
     *    *T+7.70061052263382555D-3)*T-1.71317947935716536D-2
      Y = ((((((Y*T+3.41505388391452157D-2)
     *    *T-6.04316795007737183D-2)
     *    *T+9.41616340200868389D-2)*T-1.28895621330524993D-1)
     *    *T+1.57686843969995904D-1)*T-1.86478066609466760D-1)*T +
     *    3.08508322553671039D-1
C
C     EXPANSION (0035) EVALUATED AS Y(T)  --PRECISION 19E.20
C19   Y = (((((((((((((((-1.8784179614174412800D-13
C19  *    *T+1.2089468298240983040D-12)*T-6.3074904316764487680D-12)
C19  *    *T+3.7194763684842307584D-11)*T-2.1329989100573818880D-10)
C19  *    *T+1.1532002547234216346D-9)*T-5.9434264329293442580D-9)
C19  *    *T+2.9180903864909561201D-8)*T-1.3604724152589129050D-7)
C19  *    *T+6.0059431739522308981D-7)*T-2.5029862046131434148D-6)
C19  *    *T+9.8139412868667807577D-6)*T-3.6064558781683970924D-5)
C19  *    *T+1.2368260229532060075D-4)*T-3.9393453156577578087D-4)
C19  *    *T+1.1588831957319864270D-3)*T-3.1292328666662116016D-3
C19   Y = ((((((((Y*T+7.7006105229899461374D-3)
C19  *    *T-1.7131794793558845134D-2)
C19  *    *T+3.4150538839108284692D-2)*T-6.0431679500774614881D-2)
C19  *    *T+9.4161634020088817344D-2)*T-1.2889562133052496374D-1)
C19  *    *T+1.5768684396999586206D-1)*T-1.8647806660946676075D-1)*T +
C19  *    3.0850832255367103953D-1
C
      S18AEF = G*Y
      GO TO 100
C
C     MIDDLE X
   40 G = EXP(T)
      T = 0.25D0*T - 2.0D0
C
C
C     EXPANSION (0036) EVALUATED AS Y(T)  --PRECISION 08E.09
C08   Y = (((((((((((((((-2.58491708D-6)*T+5.79902384D-6)
C08  *    *T-2.79124210D-6)*T+5.54416320D-6)*T-2.55448774D-5)
C08  *    *T+5.01438907D-5)*T-8.59465309D-5)*T+1.67961894D-4)
C08  *    *T-3.36290944D-4)*T+6.73177142D-4)*T-1.37631504D-3)
C08  *    *T+2.89353541D-3)*T-6.30122230D-3)*T+1.44861278D-2)
C08  *    *T-3.71571541D-2)*T + 1.43431782D-1
C
C     EXPANSION (0036) EVALUATED AS Y(T)  --PRECISION 09E.10
C09   Y = ((((((((((((((((+1.104771573D-6)*T-2.584917075D-6)
C09  *    *T+1.379937549D-6)*T-2.791242095D-6)*T+1.272517843D-5)
C09  *    *T-2.554487744D-5)*T+4.406764708D-5)*T-8.594653090D-5)
C09  *    *T+1.708101331D-4)*T-3.362909444D-4)*T+6.724521352D-4)
C09  *    *T-1.376315041D-3)*T+2.893626033D-3)*T-6.301222301D-3)
C09  *    *T+1.448612352D-2)*T-3.715715414D-2)*T + 1.434317819D-1
C
C     EXPANSION (0036) EVALUATED AS Y(T)  --PRECISION 12E.13
C12   Y = (((((((((((((((+2.425721534790D-8*T-6.702448677965D-8)
C12  *    *T+5.653445311965D-8)*T-1.342660594828D-7)
C12  *    *T+5.623121016284D-7)*T-1.297962120106D-6)
C12  *    *T+2.577083077721D-6)*T-5.461269051676D-6)
C12  *    *T+1.142371408329D-5)*T-2.287156222189D-5)
C12  *    *T+4.486676250931D-5)*T-8.742186248288D-5)
C12  *    *T+1.705260897922D-4)*T-3.358360330658D-4)
C12  *    *T+6.725083973161D-4)*T-1.376388812767D-3)
C12  *    *T+2.893620477892D-3
C12   Y = (((Y*T-6.301216956565D-3)
C12  *    *T+1.448612373342D-2)*T-3.715715425644D-2)*T +
C12  *    1.434317818569D-1
C
C     EXPANSION (0036) EVALUATED AS Y(T)  --PRECISION 15E.16
C15   Y = (((((((((((((((+2.825585726019810D-10
C15  *    *T-9.096828667154934D-10)*T+1.127807821205358D-9)
C15  *    *T-3.206544995575089D-9)*T+1.318013695388947D-8)
C15  *    *T-3.580576522294057D-8)*T+8.670120481253008D-8)
C15  *    *T-2.152964532328402D-7)*T+5.192357609897779D-7)
C15  *    *T-1.189131231894508D-6)*T+2.614274720462510D-6)
C15  *    *T-5.548439478014173D-6)*T+1.140334358829206D-5)
C15  *    *T-2.282804400092824D-5)*T+4.487387831119869D-5)
C15  *    *T-8.743538291189024D-5)*T+1.705245467454439D-4
C15   Y = (((((((Y*T-3.358335189754990D-4)*T+6.725085919694947D-4)
C15  *    *T-1.376389069008084D-3)*T+2.893620465323518D-3)
C15  *    *T-6.301216944612088D-3)*T+1.448612373373570D-2)
C15  *    *T-3.715715425660841D-2)*T + 1.434317818568503D-1
C
C     EXPANSION (0036) EVALUATED AS Y(T)  --PRECISION 17E.18
      Y = (((((((((((((((+2.45185252963941089D-11
     *    *T-8.46900307934754898D-11)*T+1.23188158175419302D-10)
     *    *T-3.80370174256271589D-10)*T+1.58599776268172290D-9)
     *    *T-4.66215489983794905D-9)*T+1.24131668344616429D-8)
     *    *T-3.34900221934314738D-8)*T+8.75291839187305722D-8)
     *    *T-2.17653548816447667D-7)*T+5.18632519069546106D-7)
     *    *T-1.18752840689765504D-6)*T+2.61457634142262604D-6)
     *    *T-5.54917762110482949D-6)*T+1.14032404021741277D-5)
     *    *T-2.28278155280668483D-5)*T+4.48739019580173804D-5
      Y = (((((((((Y*T-8.74354291104467762D-5)
     *    *T+1.70524543267970595D-4)
     *    *T-3.35833513200679384D-4)*T+6.72508592273773611D-4)
     *    *T-1.37638906941232170D-3)*T+2.89362046530968701D-3)
     *    *T-6.30121694459896307D-3)*T+1.44861237337359455D-2)
     *    *T-3.71571542566085323D-2)*T + 1.43431781856850311D-1
C
C     EXPANSION (0036) EVALUATED AS Y(T)  --PRECISION 19E.20
C19   Y = (((((((((((((((-4.8781177932637798400D-13
C19  *    *T+1.8591063411513098240D-12)*T-3.3263935102228889600D-12)
C19  *    *T+1.1504780908334940160D-11)*T-4.9858650699064147968D-11)
C19  *    *T+1.6385610938810420429D-10)*T-4.9724253843991245619D-10)
C19  *    *T+1.5111687324503826760D-9)*T-4.4459427230194587402D-9)
C19  *    *T+1.2503195511458724053D-8)*T-3.3744232366202447987D-8)
C19  *    *T+8.7454812402950374687D-8)*T-2.1745209965767604588D-7)
C19  *    *T+5.1867562074346417550D-7)*T-1.1876386129043432359D-6)
C19  *    *T+2.6145587489026594806D-6)*T-5.5491358755173967651D-6
C19   Y = ((((((((((((Y*T+1.1403245405046993212D-5)
C19  *    *T-2.2827826322589233732D-5)
C19  *    *T+4.4873900992550687052D-5)*T-8.7435427266706018722D-5)
C19  *    *T+1.7052454338865393184D-4)*T-3.3583351339688115207D-4)
C19  *    *T+6.7250859226473849510D-4)*T-1.3763890694005676517D-3)
C19  *    *T+2.8936204653100399418D-3)*T-6.3012169445992907185D-3)
C19  *    *T+1.4486123733735940056D-2)*T-3.7157154256608529617D-2)*T +
C19  *    1.4343178185685031071D-1
C
      S18AEF = G*Y
      GO TO 100
C
C     LARGE X
   60 G = EXP(T-LOG(T)*0.5D0)
      T = 24.0D0/T - 1.0D0
C
C
C     EXPANSION (0037) EVALUATED AS Y(T)  --PRECISION 08E.09
C08   Y = (((((+2.49574115D-8)*T+2.22642925D-7)*T+2.79704810D-6)
C08  *    *T+5.59827269D-5)*T+2.18216826D-3)*T + 4.01071065D-1
C
C     EXPANSION (0037) EVALUATED AS Y(T)  --PRECISION 09E.10
C09   Y = ((((((+3.822579830D-9)*T+2.495741155D-8)*T+2.169090557D-7)
C09  *    *T+2.797048103D-6)*T+5.598487707D-5)*T+2.182168256D-3)*T +
C09  *    4.010710651D-1
C
C     EXPANSION (0037) EVALUATED AS Y(T)  --PRECISION 12E.13
C12   Y = ((((((((((+2.340749923654D-11)*T+7.186180190509D-11)
C12  *    *T+1.579111213677D-10)*T+6.284907334603D-10)
C12  *    *T+3.440923995653D-9)*T+2.369586371205D-8)
C12  *    *T+2.171613059408D-7)*T+2.797705824730D-6)
C12  *    *T+5.598482525317D-5)*T+2.182168172172D-3)*T +
C12  *    4.010710650668D-1
C
C     EXPANSION (0037) EVALUATED AS Y(T)  --PRECISION 15E.16
C15   Y = (((((((((((((((-2.411756899398287D-13)
C15  *    *T-2.170471559606428D-12)*T-3.044730035928687D-12)
C15  *    *T+4.496440289935385D-12)*T+1.547327566378010D-11)
C15  *    *T+2.226273536199183D-11)*T+4.586796117551065D-11)
C15  *    *T+1.545697718533012D-10)*T+6.486795836711866D-10)
C15  *    *T+3.443856712516276D-9)*T+2.368833843881964D-8)
C15  *    *T+2.171604498120274D-7)*T+2.797707025776663D-6)
C15  *    *T+5.598482533619400D-5)*T+2.182168172116796D-3)*T +
C15  *    4.010710650668474D-1
C
C     EXPANSION (0037) EVALUATED AS Y(T)  --PRECISION 17E.18
      Y = (((((((((((((((-1.95679809047625728D-13
     *    *T+4.73229306831831040D-14)*T+1.44572313799118029D-12)
     *    *T+4.30812577328136192D-13)*T-4.29417106720584499D-12)
     *    *T-4.34624739357691085D-12)*T+2.82807056475555021D-12)
     *    *T+8.27719401266046976D-12)*T+1.05863621425699789D-11)
     *    *T+1.89599322920800794D-11)*T+4.82726630988879388D-11)
     *    *T+1.56147127476528831D-10)*T+6.47994117793472057D-10)
     *    *T+3.44345025431425567D-9)*T+2.36884434055843528D-8)
     *    *T+2.17160501061222148D-7)*T+2.79770701849785597D-6
      Y = ((Y*T+5.59848253337377763D-5)
     *    *T+2.18216817211694382D-3)*T +
     *    4.01071065066847416D-1
C
C     EXPANSION (0037) EVALUATED AS Y(T)  --PRECISION 19E.20
C19   Y = (((((((((((((((-3.3144842200547328000D-14
C19  *    *T+1.2008469636841472000D-14)*T+2.7186179143945420800D-13)
C19  *    *T-3.5759472348495872000D-14)*T-1.0454721503030149120D-12)
C19  *    *T-1.9782650535057817600D-13)*T+2.3978847028834467840D-12)
C19  *    *T+1.5929661492722401280D-12)*T-2.8976380738342092800D-12)
C19  *    *T-4.5692556020853309440D-12)*T-1.2796874612755333120D-12)
C19  *    *T+3.0819891159672381440D-12)*T+6.4220654560433725440D-12)
C19  *    *T+1.0449077288480794624D-11)*T+1.9688661329266574848D-11)
C19  *    *T+4.8317557581369720320D-11)*T+1.5596755827968609280D-10
C19   Y = (((((((Y*T+6.4798546599576738163D-10)
C19  *    *T+3.4434761778711018340D-9)
C19  *    *T+2.3688444310954546095D-8)*T+2.1716049913450601437D-7)
C19  *    *T+2.7977070184547602578D-6)*T+5.5984825333793869167D-5)
C19  *    *T+2.1821681721169444223D-3)*T + 4.0107106506684741582D-1
C
      S18AEF = G*Y
      GO TO 100
C
C     ERROR 1 EXIT
   80 IFAIL = P01ABF(IFAIL,1,SRNAME,0,P01REC)
      S18AEF = YBIG
C
  100 RETURN
      END
      DOUBLE PRECISION FUNCTION S18AFF(X,IFAIL)
C     MARK 7 RELEASE. NAG COPYRIGHT 1978.
C     MARK 8C REVISED. IER-269 (OCT.1980).
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C     BESSEL FUNCTION I1(X)
C
C     **************************************************************
C
C     TO EXTRACT THE CORRECT CODE FOR A PARTICULAR MACHINE-RANGE,
C     ACTIVATE THE STATEMENTS CONTAINED IN COMMENTS BEGINNING  CDD ,
C     WHERE  DD  IS THE APPROXIMATE NUMBER OF SIGNIFICANT DECIMAL
C     DIGITS REPRESENTED BY THE MACHINE
C     DELETE THE ILLEGAL DUMMY STATEMENTS OF THE FORM
C     * EXPANSION (NNNN) *
C
C     ALSO INSERT APPROPRIATE DATA STATEMENTS TO DEFINE CONSTANTS
C     WHICH DEPEND ON THE RANGE OF NUMBERS REPRESENTED BY THE
C     MACHINE, RATHER THAN THE PRECISION (SUITABLE STATEMENTS FOR
C     SOME MACHINES ARE CONTAINED IN COMMENTS BEGINNING CRD WHERE
C     D IS A DIGIT WHICH SIMPLY DISTINGUISHES A GROUP OF MACHINES).
C     DELETE THE ILLEGAL DUMMY DATA STATEMENTS WITH VALUES WRITTEN
C     *VALUE*
C
C     **************************************************************
C
C     .. Parameters ..
      CHARACTER*6                      SRNAME
      PARAMETER                        (SRNAME='S18AFF')
C     .. Scalar Arguments ..
      DOUBLE PRECISION                 X
      INTEGER                          IFAIL
C     .. Local Scalars ..
      DOUBLE PRECISION                 G, T, XBIG, XVSMAL, Y, YBIG
C     .. Local Arrays ..
      CHARACTER*1                      P01REC(1)
C     .. External Functions ..
      INTEGER                          P01ABF
      EXTERNAL                         P01ABF
C     .. Intrinsic Functions ..
      INTRINSIC                        ABS, SIGN, LOG, EXP
C     .. Data statements ..
C08   DATA XVSMAL /1.0D-4/
C09   DATA XVSMAL /3.2D-5/
C12   DATA XVSMAL /1.0D-6/
C15   DATA XVSMAL /3.2D-8/
      DATA XVSMAL /3.2D-9/
C19   DATA XVSMAL /3.2D-10/
C     FOR DEC VAX
      DATA XBIG,YBIG / 88.0D0, 1.7D+37 /
C     XBIG = LARGEST X SUCH THAT  EXP(X)/SQRT(X) .LT. MAXREAL
C             (ROUNDED DOWN)
C     YBIG = EXP(XBIG)/SQRT(XBIG)
C     FOR IEEE SINGLE PRECISION
CR0   DATA XBIG,YBIG /90.9E0,3.1E+38/
C     FOR IBM 360/370 AND SIMILAR MACHINES
CR1   DATA XBIG,YBIG /177.2D0,2.7D+75/
C     FOR DEC-10, HONEYWELL, UNIVAC 1100(S.P.)
CR2   DATA XBIG,YBIG /90.2D0,6.3D+37/
C     FOR ICL 1900
CR3   DATA XBIG,YBIG /179.3D0,2.2D+76/
C     FOR CDC 7600/CYBER
CR4   DATA XBIG,YBIG /744.9D0,4.7D+321/
C     FOR UNIVAC 1100(D.P.)
CR5   DATA XBIG,YBIG /712.3D0,3.3D+307/
C     FOR IEEE DOUBLE PRECISION
CR7   DATA XBIG,YBIG /713.0D0,1.7D+308/
C     .. Executable Statements ..
C
      T = ABS(X)
C     ERROR 1 TEST
      IF (T.GT.XBIG) GO TO 80
      IFAIL = 0
C     X RANGE TEST
      IF (T.GT.12.0D0) GO TO 60
C     TEST FOR MIDDLE RANGE
      IF (T.GT.4.0D0) GO TO 40
C     SMALL X
C     TEST FOR VERY SMALL X
      IF (T.GT.XVSMAL) GO TO 20
      S18AFF = 0.5D0*X
      GO TO 100
   20 T = 0.125D0*T*T - 1.0D0
C
C
C     EXPANSION (0043) EVALUATED AS Y(T)  --PRECISION 08E.09
C08   Y = (((((((+3.92892854D-7)*T+1.13063239D-5)*T+2.45223921D-4)
C08  *    *T+3.84762606D-3)*T+4.09286373D-2)*T+2.68657662D-1)
C08  *    *T+9.28758890D-1)*T + 1.19741655D+0
C
C     EXPANSION (0043) EVALUATED AS Y(T)  --PRECISION 09E.10
C09   Y = (((((((+3.928928537D-7)*T+1.130632386D-5)*T+2.452239209D-4)
C09  *    *T+3.847626062D-3)*T+4.092863729D-2)*T+2.686576622D-1)
C09  *    *T+9.287588901D-1)*T + 1.197416550D+0
C
C     EXPANSION (0043) EVALUATED AS Y(T)  --PRECISION 12E.13
C12   Y = (((((((((+2.330286285600D-10)*T+1.067670554920D-8)
C12  *    *T+3.923685392897D-7)*T+1.128497044797D-5)
C12  *    *T+2.452243141144D-4)*T+3.847639407499D-3)
C12  *    *T+4.092863718276D-2)*T+2.686576595217D-1)
C12  *    *T+9.287588901146D-1)*T + 1.197416549637D+0
C
C     EXPANSION (0043) EVALUATED AS Y(T)  --PRECISION 15E.16
C15   Y = ((((((((((+4.173727097882224D-12)*T+2.330286285600112D-10)
C15  *    *T+1.066627123145040D-8)*T+3.923685392897174D-7)
C15  *    *T+1.128497957799518D-5)*T+2.452243141144006D-4)
C15  *    *T+3.847639404238095D-3)*T+4.092863718276363D-2)
C15  *    *T+2.686576595220928D-1)*T+9.287588901146102D-1)*T +
C15  *    1.197416549636702D+0
C
C     EXPANSION (0043) EVALUATED AS Y(T)  --PRECISION 17E.18
      Y = (((((((((((+6.24387910353848320D-14)
     *    *T+4.17372709788222413D-12)*T+2.32856921884663846D-10)
     *    *T+1.06662712314503955D-8)*T+3.92368710996392755D-7)
     *    *T+1.12849795779951847D-5)*T+2.45224314039278904D-4)
     *    *T+3.84763940423809498D-3)*T+4.09286371827770484D-2)
     *    *T+2.68657659522092832D-1)*T+9.28758890114609554D-1)*T +
     *    1.19741654963670236D+0
C
C     EXPANSION (0043) EVALUATED AS Y(T)  --PRECISION 19E.20
C19   Y = ((((((((((((+7.9180648970649600000D-16)
C19  *    *T+6.2438791035384832000D-14)*T+4.1713516784131046400D-12)
C19  *    *T+2.3285692188466384563D-10)*T+1.0666273903797298255D-8)
C19  *    *T+3.9236871099639275503D-7)*T+1.1284979576609523316D-5)
C19  *    *T+2.4522431403927890421D-4)*T+3.8476394042384197425D-3)
C19  *    *T+4.0928637182777048352D-2)*T+2.6865765952209280388D-1)
C19  *    *T+9.2875889011460955436D-1)*T + 1.1974165496367023583D+0
C
      S18AFF = X*Y
      GO TO 100
C
C     MIDDLE X
   40 G = EXP(T)
      T = 0.25D0*T - 2.0D0
C
C
C     EXPANSION (0044) EVALUATED AS Y(T)  --PRECISION 08E.09
C08   Y = (((((((((((((((+2.04204709D-6)*T-4.39594806D-6)
C08  *    *T+1.31096677D-6)*T-1.89739179D-6)*T+1.36440625D-5)
C08  *    *T-2.21733453D-5)*T+2.35080066D-5)*T-2.41247862D-5)
C08  *    *T+5.59728936D-6)*T+8.42565753D-5)*T-3.67687996D-4)
C08  *    *T+1.17320332D-3)*T-3.40759200D-3)*T+9.76020768D-3)
C08  *    *T-2.99140925D-2)*T + 1.34142493D-1
C
C     EXPANSION (0044) EVALUATED AS Y(T)  --PRECISION 09E.10
C09   Y = ((((((((((((((((-9.013598912D-7)*T+2.042047089D-6)
C09  *    *T-7.905084934D-7)*T+1.310966772D-6)*T-7.756231083D-6)
C09  *    *T+1.364406247D-5)*T-1.721586589D-5)*T+2.350800662D-5)
C09  *    *T-2.644860471D-5)*T+5.597289356D-6)*T+8.484809269D-5)
C09  *    *T-3.676879960D-4)*T+1.173129378D-3)*T-3.407592004D-3)
C09  *    *T+9.760211205D-3)*T-2.991409248D-2)*T + 1.341424933D-1
C
C     EXPANSION (0044) EVALUATED AS Y(T)  --PRECISION 12E.13
C12   Y = (((((((((((((((-2.137968646738D-8*T+5.825193532379D-8)
C12  *    *T-4.505657546108D-8)*T+1.021824254180D-7)
C12  *    *T-4.447215247821D-7)*T+9.852042220010D-7)
C12  *    *T-1.805382790889D-6)*T+3.523606198412D-6)
C12  *    *T-6.649860183496D-6)*T+1.142085667417D-5)
C12  *    *T-1.789619803823D-5)*T+2.473724235346D-5)
C12  *    *T-2.620655968490D-5)*T+5.217818255143D-6)
C12  *    *T+8.480011929240D-5)*T-3.676264096858D-4)
C12  *    *T+1.173134117233D-3
C12   Y = (((Y*T-3.407596468626D-3)
C12  *    *T+9.760211025571D-3)*T-2.991409238989D-2)*T +
C12  *    1.341424932927D-1
C
C     EXPANSION (0044) EVALUATED AS Y(T)  --PRECISION 15E.16
C15   Y = (((((((((((((((-2.585602855843796D-10
C15  *    *T+8.263393060071838D-10)*T-9.914352278371534D-10)
C15  *    *T+2.772708966945032D-9)*T-1.146662778794534D-8)
C15  *    *T+3.062872297108590D-8)*T-7.213105372475229D-8)
C15  *    *T+1.741345209187592D-7)*T-4.060152623900778D-7)
C15  *    *T+8.884275924436715D-7)*T-1.838822740857724D-6)
C15  *    *T+3.601181743921443D-6)*T-6.631537188516544D-6)
C15  *    *T+1.138210964994464D-5)*T-1.790260037692807D-5)
C15  *    *T+2.474928446808203D-5)*T-2.620517107156340D-5
C15   Y = (((((((Y*T+5.215578508861103D-6)*T+8.479994409371570D-5)
C15  *    *T-3.676261813644457D-4)*T+1.173134128546845D-3)
C15  *    *T-3.407596479277477D-3)*T+9.760211025286696D-3)
C15  *    *T-2.991409238974067D-2)*T + 1.341424932926982D-1
C
C     EXPANSION (0044) EVALUATED AS Y(T)  --PRECISION 17E.18
      Y = (((((((((((((((-2.27061376122617856D-11
     *    *T+7.79929176497056645D-11)*T-1.10970391104678003D-10)
     *    *T+3.38883570696523350D-10)*T-1.41575617446629553D-9)
     *    *T+4.11321223904934809D-9)*T-1.07563514207617768D-8)
     *    *T+2.84961041291017650D-8)*T-7.28978293484163628D-8)
     *    *T+1.76305222240064495D-7)*T-4.05456611578551130D-7)
     *    *T+8.86951515545183908D-7)*T-1.83910206626348772D-6)
     *    *T+3.60186151617732531D-6)*T-6.63144162982509821D-6)
     *    *T+1.13818992442463952D-5)*T-1.79026222757948636D-5
      Y = (((((((((Y*T+2.47493270133518925D-5)
     *    *T-2.62051678511418163D-5)
     *    *T+5.21557319070236939D-6)*T+8.47999438119288094D-5)
     *    *T-3.67626180992174570D-4)*T+1.17313412855965374D-3)
     *    *T-3.40759647928956354D-3)*T+9.76021102528646704D-3)
     *    *T-2.99140923897405570D-2)*T + 1.34142493292698178D-1
C
C     EXPANSION (0044) EVALUATED AS Y(T)  --PRECISION 19E.20
C19   Y = (((((((((((((((+4.5775493193216819200D-13
C19  *    *T-1.7377647145171025920D-12)*T+3.0683838993620008960D-12)
C19  *    *T-1.0541784610642067456D-11)*T+4.5665794931233390592D-11)
C19  *    *T-1.4898399423473962189D-10)*T+4.4748047024971893965D-10)
C19  *    *T-1.3458111447069821501D-9)*T+3.9122095980161710162D-9)
C19  *    *T-1.0840504034690950758D-8)*T+2.8732496785343932531D-8)
C19  *    *T-7.2828311971692262588D-8)*T+1.7611785989651199931D-7)
C19  *    *T-4.0549690005824350650D-7)*T+8.8705402749571279009D-7)
C19  *    *T-1.8390856219860622587D-6)*T+3.6018226816369370946D-6
C19   Y = ((((((((((((Y*T-6.6314463061664910712D-6)
C19  *    *T+1.1381909286711502521D-5)
C19  *    *T-1.7902621373343015831D-5)*T+2.4749325297973732372D-5)
C19  *    *T-2.6205167963948297248D-5)*T+5.2155733732522876693D-6)
C19  *    *T+8.4799943820374214419D-5)*T-3.6762618100311115229D-4)
C19  *    *T+1.1731341285593238429D-3)*T-3.4075964792892586671D-3)
C19  *    *T+9.7602110252864721110D-3)*T-2.9914092389740559536D-2)*T +
C19  *    1.3414249329269817831D-1
C
      S18AFF = SIGN(G*Y,X)
      GO TO 100
C
C     LARGE X
   60 G = EXP(T-LOG(T)*0.5D0)
      T = 24.0D0/T - 1.0D0
C
C
C     EXPANSION (0045) EVALUATED AS Y(T)  --PRECISION 08E.09
C08   Y = (((((-2.96738127D-8)*T-2.78984102D-7)*T-3.82721490D-6)
C08  *    *T-9.12451446D-5)*T-6.40545370D-3)*T + 3.92624494D-1
C
C     EXPANSION (0045) EVALUATED AS Y(T)  --PRECISION 09E.10
C09   Y = ((((((-4.383292624D-9)*T-2.967381267D-8)*T-2.724091627D-7)
C09  *    *T-3.827214904D-6)*T-9.124761019D-5)*T-6.405453697D-3)*T +
C09  *    3.926244942D-1
C
C     EXPANSION (0045) EVALUATED AS Y(T)  --PRECISION 12E.13
C12   Y = ((((((((((-2.533926898209D-11)*T-7.774606791952D-11)
C12  *    *T-1.734215243288D-10)*T-7.068696867564D-10)
C12  *    *T-3.965182881278D-9)*T-2.826186206848D-8)
C12  *    *T-2.726853285086D-7)*T-3.827950034494D-6)
C12  *    *T-9.124755347160D-5)*T-6.405453603544D-3)*T +
C12  *    3.926244942041D-1
C
C     EXPANSION (0045) EVALUATED AS Y(T)  --PRECISION 15E.16
C15   Y = (((((((((((((((+2.864463574710190D-13)
C15  *    *T+2.296326081181966D-12)*T+3.031081657644896D-12)
C15  *    *T-4.930682983830821D-12)*T-1.625784942705443D-11)
C15  *    *T-2.360757461732030D-11)*T-4.984995964878307D-11)
C15  *    *T-1.704720475191773D-10)*T-7.288220600471799D-10)
C15  *    *T-3.967981996319198D-9)*T-2.825360264463809D-8)
C15  *    *T-2.726844939064504D-7)*T-3.827951362129912D-6)
C15  *    *T-9.124755355333729D-5)*T-6.405453603482220D-3)*T +
C15  *    3.926244942041166D-1
C
C     EXPANSION (0045) EVALUATED AS Y(T)  --PRECISION 17E.18
      Y = (((((((((((((((+1.99448557598015488D-13
     *    *T-5.77176811730370560D-14)*T-1.48765082315961139D-12)
     *    *T-3.95353303949377536D-13)*T+4.47735589657057690D-12)
     *    *T+4.42966462319664333D-12)*T-3.05957293450420224D-12)
     *    *T-8.69631766630563635D-12)*T-1.11795516742222899D-11)
     *    *T-2.02947854602758139D-11)*T-5.23524129533553498D-11)
     *    *T-1.72060490748583241D-10)*T-7.28107961041827952D-10)
     *    *T-3.96757162863209348D-9)*T-2.82537120880041703D-8)
     *    *T-2.72684545741400871D-7)*T-3.82795135453556215D-6
      Y = ((Y*T-9.12475535508497109D-5)
     *    *T-6.40545360348237412D-3)*T +
     *    3.92624494204116555D-1
C
C     EXPANSION (0045) EVALUATED AS Y(T)  --PRECISION 19E.20
C19   Y = (((((((((((((((+3.3679618863005696000D-14
C19  *    *T-1.3569887628689408000D-14)*T-2.7733286267020902400D-13)
C19  *    *T+4.6135638859186176000D-14)*T+1.0726764420299489280D-12)
C19  *    *T+1.7181030171593932800D-13)*T-2.4822246330641940480D-12)
C19  *    *T-1.5863951335214284800D-12)*T+3.0522499336199864320D-12)
C19  *    *T+4.7030271230869504000D-12)*T+1.2431528563850690560D-12)
C19  *    *T-3.2807424254310154240D-12)*T-6.7639380417875619840D-12)
C19  *    *T-1.1056476687880072192D-11)*T-2.1055276286591258112D-11)
C19  *    *T-5.2393348632806530560D-11)*T-1.7187282226057648858D-10
C19   Y = (((((((Y*T-7.2809998359709274272D-10)
C19  *    *T-3.9675987527407277127D-9)
C19  *    *T-2.8253712929434771752D-8)*T-2.7268454372357710872D-7)
C19  *    *T-3.8279513544952731211D-6)*T-9.1247553550908501103D-5)
C19  *    *T-6.4054536034823746802D-3)*T + 3.9262449420411655531D-1
C
      S18AFF = SIGN(G*Y,X)
      GO TO 100
C
C     ERROR 1 EXIT
   80 IFAIL = P01ABF(IFAIL,1,SRNAME,0,P01REC)
      S18AFF = YBIG
C
  100 RETURN
      END
      REAL FUNCTION TG01BD(IX,N,U,S,D,X)
C###### 04/08/70LAST LIBRARY UPDATE
C
C**********************************************************************
C
C      TG01BD - FUNCTION ROUTINE TO EVALUATE A CUBIC SPLINE GIVEN SPLINE
C     VALUES AND FIRST DERIVATIVE VALUES AT THE GIVEN KNOTS.
C
C     THE SPLINE VALUE IS DEFINED AS ZERO OUTSIDE THE KNOT RANGE,WHICH
C     IS EXTENDED BY A ROUNDING ERROR FOR THE PURPOSE.
C
C                 F = TG01BD(IX,N,U,S,D,X)
C
C       IX    ALLOWS CALLER TO TAKE ADVANTAGE OF SPLINE PARAMETERS SET
C             ON A PREVIOUS CALL IN CASES WHEN X POINT FOLLOWS PREVIOUS
C             X POINT. IF IX < 0 THE WHOLE RANGE IS SEARCHED FOR KNOT
C             INTERVAL; IF IX > 0 IT IS ASSUMED THAT X IS GREATER THAN
C             THE X OF THE PREVIOUS CALL AND SEARCH STARTED FROM THERE.
C       N     NUMBER OF KNOTS.
C       U     THE KNOTS.
C       S     THE SPLINE VALUES.
C       D     THE FIRST DERIVATIVE VALUES OF THE SPLINE AT THE KNOTS.
C       X     THE POINT AT WHICH THE SPLINE VALUE IS REQUIRED.
C       F     THE VALUE OF THE SPLINE AT THE POINT X.
C
C                                      MODIFIED JULY 1970
C
C**********************************************************************
C
      IMPLICIT REAL*8(A-H,O-Z)
C
C     ALLOWABLE ROUNDING ERROR ON POINTS AT EXTREAMS OF KNOT RANGE
C     IS 2**IEPS*MAX(!U(1)!,!U(N)!).
      INTEGER*4 IFLG/0/,IEPS/-50/
      DIMENSION U(1),S(1),D(1)
C
C       TEST WETHER POINT IN RANGE.
      IF(X.LT.U(1)) GO TO 990
      IF(X.GT.U(N)) GO TO 991
C
C       JUMP IF KNOT INTERVAL REQUIRES RANDOM SEARCH.
      IF(IX.LT.0.OR.IFLG.EQ.0) GO TO 12
C       JUMP IF KNOT INTERVAL SAME AS LAST TIME.
      IF(X.LE.U(J+1)) GO TO 8
C       LOOP TILL INTERVAL FOUND.
    1 J=J+1
   11 IF(X.GT.U(J+1)) GO TO 1
      GO TO 7
C
C       ESTIMATE KNOT INTERVAL BY ASSUMING EQUALLY SPACED KNOTS.
   12 J=DABS(X-U(1))/(U(N)-U(1))*(N-1)+1
C       ENSURE CASE X=U(N) GIVES J=N-1.
      J=MIN0(J,N-1)
C       INDICATE THAT KNOT INTERVAL INSIDE RANGE HAS BEEN USED.
      IFLG=1
C       SEARCH FOR KNOT INTERVAL CONTAINING X.
      IF(X.GE.U(J)) GO TO 11
    2 J=J-1
      IF(X.LT.U(J)) GO TO 2
C
C       CALCULATE SPLINE PARAMETERS FOR JTH INTERVAL.
    7 H=U(J+1)-U(J)
      Q1=H*D(J)
      Q2=H*D(J+1)
      SS=S(J+1)-S(J)
      B=3D0*SS-2D0*Q1-Q2
      A=Q1+Q2-2D0*SS
C
C       CALCULATE SPLINE VALUE.
    8 Z=(X-U(J))/H
      TG01BD=((A*Z+B)*Z+Q1)*Z+S(J)
      RETURN
C       TEST IF X WITHIN ROUNDING ERROR OF U(1).
  990 IF(X.LE.U(1)-2D0**IEPS*DMAX1(DABS(U(1)),DABS(U(N)))) GO TO 99
      J=1
      GO TO 7
C       TEST IF X WITHIN ROUNDING ERROR OF U(N).
  991 IF(X.GE.U(N)+2D0**IEPS*DMAX1(DABS(U(1)),DABS(U(N)))) GO TO 99
      J=N-1
      GO TO 7
   99 IFLG=0
C       FUNCTION VALUE SET TO ZERO FOR POINTS OUTSIDE THE RANGE.
      TG01BD=0D0
      RETURN
      END
      SUBROUTINE VB06AD(M,N,XD,YD,WD,RD,XN,FN,GN,DN,THETA,IPRINT,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      LOGICAL SWITCH
      DIMENSION XD(1),YD(1),WD(1),RD(1),XN(1),FN(1),GN(1),DN(1),
     1THETA(1),W(1)
C     SET THE RELATIVE ACCURACY OF THE COMPUTER ARITHMETIC
      PREC=16.0D0**(-14)
C     RESERVE THE FIRST FIFTEEN LOCATIONS OF W FOR
C     COMPUTED VALUES AND THIRD DERIVATIVES OF B-SPLINES
      W(1)=0.D0
      W(7)=0.D0
C     THEN SET THE KNOT POSITIONS INCLUDING SIX OUTSIDE THE RANGE
C     AND DELETE ANY KNOTS NOT IN STRICTLY ASCENDING ORDER
      I=1
      NN=1
   10 W(NN+18)=XN(I)
   20 IF (I.GE.N) GO TO 30
      I=I+1
      IF (XN(I).LE.W(NN+18)) GO TO 20
      NN=NN+1
      XN(NN)=XN(I)
      THETA(NN)=THETA(I)
      GO TO 10
   30 IF (NN.GE.N) GO TO 50
      N=NN
      PRINT 40
   40 FORMAT (//5X,'SOME KNOTS HAVE BEEN DELETED BY VB06AD')
   50 IF (NN.GE.2) GO TO 70
      PRINT 60
   60 FORMAT (//5X,'THE EXECUTION OF VB06AD HAS BEEN STOPPED ',
     1'BECAUSE THERE ARE TOO FEW KNOTS')
      GO TO 770
   70 DO 80 I=1,3
      W(19-I)=2.D0*W(20-I)-W(21-I)
   80 W(NN+I+18)=2.D0*W(NN+I+17)-W(NN+I+16)
C     DELETE ANY DATA POINTS THAT ARE NOT IN ASCENDING ORDER
      I=1
      MM=1
      XD(MM)=DMIN1(DMAX1(XD(1),XN(1)),XN(NN))
   90 IF (I.GE.M) GO TO 100
      I=I+1
      IF (XD(I).LT.XD(MM)) GO TO 90
      MM=MM+1
      XD(MM)=DMIN1(XD(I),XN(NN))
      YD(MM)=YD(I)
      WD(MM)=WD(I)
      GO TO 90
  100 IF (MM.GE.M) GO TO 120
      M=MM
      PRINT 110
  110 FORMAT (//5X,'SOME DATA POINTS HAVE BEEN DELETED BY VB06AD')
C     COMPUTE THE THIRD DERIVATIVES OF THE B-SPLINES AT XN(1)
  120 SWITCH=.TRUE.
      IN=1
      GO TO 280
  130 DO 140 I=2,5
      W(I+6)=W(I)
  140 W(I+10)=W(I)
C     INITIALIZE THE FACTORIZATION OF THE LEAST SQUARES MATRIX
      KS=NN+22
      KU=KS+22
      DO 150 I=KS,KU
  150 W(I)=0.D0
      IM=1
      IN=2
      KK=KS
      SWITCH=.FALSE.
C     TEST WHETHER A DATA POINT OR A KNOT DEFINES THE NEXT EQUATION
  160 KR=KK
      NC=5
      IF (IM.GT.MM) GO TO 260
      IF (XD(IM).GT.XN(IN)) GO TO 260
C     CALCULATE THE VALUES OF THE B-SPLINES AT XD(IM)
  170 W(2)=1.D0/(W(IN+18)-W(IN+17))
      DO 190 J=2,4
      I=J+1
      W(I)=0.D0
  180 II=IN+I+16
      C=(XD(IM)-W(II-J))*W(I-1)+(W(II)-XD(IM))*W(I)
      IF (DABS(C).LE.PREC) C=0.
      W(I)=C/(W(II)-W(II-J))
      I=I-1
      IF (I.GE.2) GO TO 180
  190 CONTINUE
      IF (SWITCH) GO TO 540
C     COMPLETE THE EQUATION OF THE DATA POINT
      W(6)=YD(IM)
      C=WD(IM)**2
      IM=IM+1
C     REVISE NC PRIOR TO THE NEXT GIVENS TRANSFORMATION
  200 NC=NC-1
      IF (NC.LE.0) GO TO 160
C     MAKE THE GIVENS TRANSFORMATION
  210 ALPHA=C*W(2)
      IF (ALPHA.NE.0.D0) GO TO 230
      DO 220 I=1,NC
  220 W(I+1)=W(I+2)
      GO TO 250
  230 DSTAR=W(KR)+ALPHA*W(2)
      BETA=W(KR)/DSTAR
      C=BETA*C
      GAMMA=ALPHA/DSTAR
      W(KR)=DSTAR
      DSTAR=W(2)
      DO 240 I=1,NC
      W(I+1)=W(I+2)-DSTAR*W(KR+I)
  240 W(KR+I)=BETA*W(KR+I)+GAMMA*W(I+2)
  250 KR=KR+7
      GO TO 200
C     ADVANCE THE CURRENT ROW OF THE UPPER TRIANGULAR MATRIX
  260 IF (IN.GE.NN) GO TO 330
      W(KK+28)=0.D0
      DO 270 I=4,28,6
      W(KK+I+1)=W(KK+I)
  270 W(KK+I)=0.D0
      KK=KK+7
C     CALCULATE THE THIRD DERIVATIVES OF THE B-SPLINES AT XN(IN+)
  280 W(2)=6.D0/(W(IN+19)-W(IN+18))
      DO 300 J=2,4
      I=J+1
      W(I)=0.D0
  290 II=IN+I+17
      W(I)=(W(I-1)-W(I))/(W(II)-W(II-J))
      I=I-1
      IF (I.GE.2) GO TO 290
  300 CONTINUE
      IF (SWITCH) GO TO 130
C     SET UP THE EQUATION FOR THE SMOOTHING TERM AT A KNOT
      ALPHA=0.D0
      I=6
  310 W(I)=W(I-1)-ALPHA
      I=I-1
      IF (I.LE.1) GO TO 320
      ALPHA=W(I+6)
      W(I+6)=W(I)
      GO TO 310
  320 C=THETA(IN)**2
      IN=IN+1
      GO TO 210
C     BACK-SUBSTITUTE TO OBTAIN THE MULTIPLIERS OF THE B-SPLINES
  330 KU=KK
      KUU=KU+21
      KB=KS+7
      IBS=5
      DO 335 I=3,6
      KK=KK+7
  335 W(KK-2)=W(KK-I)
  340 KK=KUU+IBS
      KR=KK
  350 KK=KK-7
      K=KK
      J=KK-IBS
  360 J=J+1
      K=K+7
      W(KK)=W(KK)-W(J)*W(K)
      IF (K.LT.KR) GO TO 360
      KR=MIN0(KR,KK+21)
      IF (KK.GT.KB) GO TO 350
      IF (IBS-6) 370,470,500
  370 SWITCH=.TRUE.
      ALPHA=1.D0
      DO 380 K=KS,KUU,7
  380 ALPHA=DMIN1(ALPHA,W(K))
      IF (ALPHA.LE.0.D0) GO TO 592
      IF (NN.LE.2) GO TO 600
C     NEXT A SPLINE WITH FIXED VALUES OF S'''(XN(1)) AND S'''(XN(N))
C     IS CALCULATED TO OBTAIN THE OPTIMAL EXTREME VALUES OF S'''(X).
C     SET THE COEFFICIENTS OF S'''(X(1)) AND S'''(X(N))
      KR=KS
      KK=KU
      DO 410 I=8,11
      W(KK+4)=W(I)
      KK=KK+7
      W(KR+6)=W(I+4)
  410 KR=KR+7
      DO 420 K=KR,KUU,7
  420 W(K+6)=0.D0
      IBS=6
      KRR=KS
C     BACK-SUBSTITUTE USING THE TRANSPOSE OF THE UPPER TRIANGULAR
C     MATRIX TO CALCULATE THE RESIDUALS OF EACH EXTRA SPLINE
  425 I=KRR+7
      KR=KRR
      DO 450 KK=I,KUU,7
      J=KK
      K=KK
  430 K=K-7
      J=J-6
      W(KK+IBS)=W(KK+IBS)-W(J)*W(K+IBS)
      IF (K.GT.KR) GO TO 430
  450 KR=MAX0(KR,KK-21)
      DO 460 KK=KRR,KUU,7
  460 W(KK+IBS)=W(KK+IBS)/W(KK)
      IF (IBS.EQ.6) GO TO 340
      DO 465 K=KS,KUU,7
      W(K+7)=0.D0
  465 IF (K.GE.KU) W(K+7)=W(K+4)
      IBS=7
      GO TO 340
C     BACK-SUBSTITUTE TO GIVE THE VECTOR FOR S'''(X(N))
  470 IBS=4
      KRR=KU
      GO TO 425
C     AT EACH DATA POINT CALCULATE THE RESIDUALS OF THREE SPLINES
  500 DO 520 I=6,10
  520 W(I)=0.D0
      IN=2
      KK=KS
      DO 580 IM=1,MM
  530 IF (XD(IM).LE.XN(IN)) GO TO 170
      IN=IN+1
      KK=KK+7
      GO TO 530
  540 KR=KK
      W(11)=YD(IM)
      W(12)=0.D0
      W(13)=0.D0
      DO 560 I=2,5
      DO 550 J=5,7
  550 W(J+6)=W(J+6)-W(I)*W(KR+J)
  560 KR=KR+7
      DO 570 I=11,13
  570 W(I)=WD(IM)*W(I)
C     FORM THE NORMAL EQUATIONS TO FIX S'''(X) AT THE ENDS OF THE RANGE
      DO 580 I=12,13
      K=I+I-29
      DO 580 J=11,I
  580 W(K+J)=W(K+J)+W(I)*W(J)
C     CALCULATE THE B-SPLINE MULTIPLIERS OF THE REQUIRED FIT
      ALPHA=(W(8)*W(9)-W(6)*W(10))/(W(7)*W(10)-W(9)**2)
      BETA=(-W(8)-ALPHA*W(9))/W(10)
      DO 590 K=KS,KUU,7
  590 W(K+5)=W(K+5)+ALPHA*W(K+6)+BETA*W(K+7)
      GO TO 600
C     PRINT A DIAGNOSTIC MESSAGE IF THERE IS INSUFFICIENT DATA
  592 PRINT 594
  594 FORMAT (//,'::',5X,'THE RESULTS FROM VB06AD ARE ',
     .'UNRELIABLE BECAUSE ',
     1'THERE IS INSUFFICIENT DATA TO DEFINE THE SPLINE UNIQUELY')
C     CALCULATE THE REQUIRED VALUES OF S(X) AND S'(X) AT THE KNOTS
  600 KK=KS
      DO 620 IN=1,NN
      FN(IN)=0.D0
      GN(IN)=0.D0
      W(2)=(W(IN+19)-W(IN+18))/((W(IN+19)-W(IN+17))*(W(IN+19)-W(IN+16)))
      W(3)=(W(IN+18)-W(IN+17))/((W(IN+19)-W(IN+17))*(W(IN+20)-W(IN+17)))
      W(4)=0.D0
      DO 610 I=1,3
      J=IN+I+14
      ALPHA=((XN(IN)-W(J))*W(I)+(W(J+4)-XN(IN))*W(I+1))/(W(J+4)-W(J))
      FN(IN)=FN(IN)+ALPHA*W(KK+5)
      BETA=3.D0*(W(I)-W(I+1))/(W(J+4)-W(J))
      GN(IN)=GN(IN)+BETA*W(KK+5)
  610 KK=KK+7
  620 KK=KK-14
C     CALCULATE THE THIRD DERIVATIVE DISCONTINUITIES OF THE SPLINE
      DN(1)=0.D0
      IN=1
      IM=1
  630 H=XN(IN+1)-XN(IN)
      ALPHA=FN(IN)-FN(IN+1)+H*GN(IN)
      BETA=ALPHA+FN(IN)-FN(IN+1)+H*GN(IN+1)
      DN(IN+1)=6.D0*BETA/H**3
      DN(IN)=DN(IN+1)-DN(IN)
      IN=IN+1
C     CALCULATE THE RESIDUALS AT THE DATA POINTS
  640 IF (IM.GT.MM) GO TO 650
      IF (XD(IM).GT.XN(IN)) GO TO 630
      C=(XD(IM)-XN(IN-1))/H
      RD(IM)=YD(IM)-C*FN(IN)+(C-1.D0)*(FN(IN-1)+C*(ALPHA-C*BETA))
      IM=IM+1
      GO TO 640
  650 IF (IN.LT.NN) GO TO 630
C     PROVIDE PRINTING IF REQUESTED
      IF (IPRINT.LE.0) GO TO 770
      PRINT 660
  660 FORMAT (1H1,35X,'SPLINE APPROXIMATION OBTAINED BY VB06AD'//4X,'I',
     19X,'XN(I)',18X,'FN(I)',18X,'GN(I)',13X,'3RD DERIV CHANGE',
     211X,'THETA(I)'//)
      I=1
  670 PRINT 680,I,XN(I),FN(I),GN(I)
  680 FORMAT (I5,5D23.14)
  690 I=I+1
      IF (I-NN) 700,670,710
  700 PRINT 680,I,XN(I),FN(I),GN(I),DN(I),THETA(I)
      GO TO 690
  710 PRINT 720
  720 FORMAT (///4X,'I',9X,'XD(I)',18X,'YD(I)',18X,'WD(I)',19X,'FIT',
     118X,'RESIDUAL'//)
      IN=2
      DO 760 IM=1,MM
  730 IF (XD(IM).LE.XN(IN)) GO TO 750
      IN=IN+1
      IF (SWITCH) GO TO 730
      SWITCH=.TRUE.
      PRINT 740
  740 FORMAT (5X)
      GO TO 730
  750 C=YD(IM)-RD(IM)
      PRINT 680,IM,XD(IM),YD(IM),WD(IM),C,RD(IM)
  760 SWITCH=.FALSE.
  770 RETURN
      END
      SUBROUTINE VC03AD (M,N,XD,YD,WD,RD,XN,FN,GN,DN,THETA,IPRINT,W)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION XD(1),YD(1),WD(1),RD(1),XN(1),FN(1),GN(1),DN(1),
     1THETA(1),W(1)
      NDIMAX=N
C     OMIT DATA WITH ZERO WEIGHTS
      MM=1
      J=1
      HLF=0.5D0
      DO 1 I=2,M
      IF (WD(I)) 3,2,3
    2 IF (I-M) 4,3,3
    4 W(J)=XD(I)
      W(J+1)=YD(I)
      W(J+2)=DFLOAT(I)
      J=J+3
      GO TO 1
    3 MM=MM+1
      XD(MM)=XD(I)
      YD(MM)=YD(I)
      WD(MM)=WD(I)
    1 CONTINUE
      J=1
      K=MM
    7 IF (K-M) 5,6,6
    5 K=K+1
      XD(K)=W(J)
      YD(K)=W(J+1)
      WD(K)=W(J+2)
      J=J+3
      GO TO 7
C     INITIALIZATION OF ITERATIONS
    6 JA=1
      IF (WD(1)) 9,8,9
    8 JA=2
    9 N=5
      IF(N-NDIMAX)100,100,110
  100 XN(1)=XD(1)
      XN(5)=XD(MM)
      XN(3)=HLF*(XN(1)+XN(5))
      XN(2)=HLF*(XN(1)+XN(3))
      XN(4)=HLF*(XN(3)+XN(5))
      IW=7*MM+46
      DO 10 I=1,5
      J=IW+I
      W(J)=XN(I)
   10 CONTINUE
      IP=-IPRINT
C     CALCULATE THE HISTOGRAMS FOR NEW SCALE FACTORS
   11 J=0
      K=1
      SA=0.0D0
   12 SA=SA+WD(K)**2
      K=K+1
      IF (XD(K)-XD(1)) 12,12,13
   13 SA=(SA+HLF*WD(K)**2)/(XD(K)-XD(1))
   14 J=J+1
      IF (XD(K)-XN(J+1)) 15,15,16
   16 GN(J)=SA*(XN(J+1)-XN(J))
      GO TO 14
   15 GN(J)=SA*(XD(K)-XN(J))+HLF*WD(K)**2
   17 K=K+1
      IF (K-MM) 18,18,19
   18 IF (XD(K)-XN(J+1)) 20,20,21
   20 GN(J)=GN(J)+WD(K)**2
      GO TO 17
   21 SA=HLF*(WD(K-1)**2+WD(K)**2)/(XD(K)-XD(K-1))
      GN(J)=GN(J)-HLF*WD(K-1)**2+SA*(XN(J+1)-XD(K-1))
      GO TO 14
C     CALCULATE THE NEW SCALE FACTORS
   19 K=IW+2
      GN(1)=0.00025216D0*GN(1)*(W(K)-W(K-1))**8/(XN(2)-XN(1))
      DO 22 J=3,N
      IF (XN(J)-W(K)) 23,23,24
   24 K=K+1
   23 GN(J-1)=0.00025216D0*GN(J-1)*(W(K)-W(K-1))**8/(XN(J)-XN(J-1))
      GN(J-2)=DLOG(GN(J-2)+GN(J-1))
   22 CONTINUE
      NN=N-2
      HS=1.386294D0
      DO 97 J=2,NN
      GN(J)=DMIN1(GN(J),GN(J-1)+HS)
   97 CONTINUE
      J=NN-1
   98 GN(J)=DMIN1(GN(J),GN(J+1)+HS)
      J=J-1
      IF (J) 99,99,98
   99 DO 25 J=3,N
      THETA(J-1)=DSQRT(DEXP(GN(J-2))/(XN(J)-XN(J-2)))
   25 CONTINUE
C     CALCULATE THE SPLINE APPROXIMATION WITH CURRENT KNOTS
      CALL VB06AD (MM,N,XD,YD,WD,RD,XN,FN,GN,DN,THETA,IP,W)
C     APPLY STATISTICAL TEST FOR EXTRA KNOTS
      J=IW+1
      JJ=0
      K=JA
      TMAX=0.0D0
      IIS=1
   26 KC=-1
      SW=0.0D0
      SR=0.0D0
      RP=0.0D0
      J=J+1
      JJ=JJ+1
      W(JJ)=0.0D0
   27 IF (W(J)-XD(K)) 28,29,30
   30 KC=KC+1
      SW=SW+RD(K)**2
      SR=SR+RP*RD(K)
      RP=RD(K)
      K=K+1
      GO TO 27
   29 IF (WD(K)) 31,28,31
   31 KC=KC+1
      SW=SW+RD(K)**2
      SR=SR+RP*RD(K)
   28 IF (SR) 32,32,33
   33 SW=(SW/SR)**2
      RP=DSQRT(SR/DFLOAT(KC))
      IF (DFLOAT(KC)-SW) 32,32,34
   34 GO TO (35,36,37),IIS
   35 PRP=RP
      IIS=3
      IF (DFLOAT(KC)-2.D0*SW) 38,38,39
   37 W(JJ-1)=PRP
      TMAX=DMAX1(TMAX,PRP)
   39 IIS=2
   36 W(JJ)=RP
      TMAX=DMAX1(TMAX,RP)
      GO TO 38
   32 IIS=1
   38 IF (W(J)-XN(N)) 26,40,40
C     TEST WHETHER ANOTHER ITERATION IS REQUIRED
   40 IF (TMAX) 41,41,42
C     CALCULATE NEW TREND ARRAY, INCLUDING LARGER TRENDS ONLY
   42 TMAX=HLF*TMAX
      I=0
      J=1
      JW=1
      K=IW+1
      THETA(JW)=W(K)
   43 I=I+1
      K=K+1
      IF (W(I)-TMAX) 44,44,45
   44 JW=JW+1
      THETA(JW)=W(K)
   46 FN(J)=0.0D0
      J=J+1
      IF (W(K)-XN(J)) 47,47,46
   45 JW=JW+2
      THETA(JW-1)=HLF*(W(K-1)+W(K))
      THETA(JW)=W(K)
      IF (XN(J+1)-THETA(JW-1)) 46,46,48
   48 FN(J)=1.0D0
      J=J+1
   47 IF (J-N) 43,49,49
C     MAKE KNOT SPACINGS BE USED FOUR TIMES
   49 IK=1
      KL=1
      FN(2)=DMAX1(FN(1),FN(2))
      GO TO 102
   50 K=KL+3
   51 IF (FN(K)) 52,52,53
   52 K=K-1
      IF (K-KL) 74,74,51
   53 K=K-1
      FN(K)=1.0D0
      IF (K-KL) 74,74,53
  102 K=KL+3
   54 K=K+1
      IF (K-N) 55,56,56
   55 IF (XN(K+1)-XN(K)-1.5*(XN(K)-XN(K-1))) 54,54,56
   56 KU=K
      FN(K-2)=DMAX1(FN(K-2),FN(K-1))
   57 KKU=K
   58 K=K-1
      IF (K-KL) 59,59,60
   60 IF (XN(K)-XN(K-1)-1.5*(XN(K+1)-XN(K))) 58,58,61
   61 FN(K+1)=DMAX1(FN(K),FN(K+1))
   59 KKL=K
      KZ=4
      IF (FN(K)) 62,62,63
   63 K=K+1
      IF (K-KKU) 64,65,65
   64 IF (FN(K)) 66,66,63
   66 KZ=0
   62 KZ=KZ+1
      K=K+1
      IF (K-KKU) 67,65,65
   67 IF (FN(K)) 62,62,68
   68 IF (KZ-3) 69,69,70
   69 J=K-KZ
   71 FN(J)=1.0D0
      J=J+1
      IF (J-K) 71,63,63
   70 IF (K+1-KKU) 72,65,65
   72 K=K+1
      FN(K)=1.0D0
      GO TO 63
   65 IF (KL-KKL) 73,50,50
   73 FN(KKL-2)=DMAX1(FN(KKL-2),FN(KKL+1))
      FN(KKL-1)=DMAX1(FN(KKL-1),FN(KKL+3))
   75 K=KKL-4
   78 IF (FN(K)) 76,76,77
   76 K=K+1
      IF (K-KKL) 78,79,79
   77 FN(K)=1.0D0
      K=K+1
      IF (K-KKL) 77,79,79
   79 GO TO (57,80),IK
   74 IF (KU-N) 81,82,82
   81 KL=KU
      FN(KL+1)=DMAX1(FN(KL+1),FN(KL-2))
      FN(KL)=DMAX1(FN(KL),FN(KL-4))
      GO TO 102
   82 IK=2
      KKL=N
      GO TO 75
C     INSERT EXTRA KNOTS FOR NEW APPROXIMATION
   80 DO 83 J=1,N
      GN(J)=XN(J)
   83 CONTINUE
      NN=1
      DO 84 J=2,N
      IF (FN(J-1)) 85,85,86
   86 NN=NN+1
      XN(NN)=HLF*(GN(J-1)+GN(J))
   85 NN=NN+1
      XN(NN)=GN(J)
   84 CONTINUE
      IF(N-NDIMAX)101,110,110
  110 PRINT 111,N
  111 FORMAT(///' ARRAY SIZES TOO SMALL.  N =',I6,///)
      N=-N
      GO TO 90
  101 N=NN
      IW=7*MM+8*N+6
      DO 87 J=1,JW
      I=IW+J
      W(I)=THETA(J)
   87 CONTINUE
      GO TO 11
C     RESTORE DATA WITH ZERO WEIGHTS
   41 IF (MM-M) 88,89,89
   89 IF (IPRINT) 90,90,91
   88 J=-2
      K=MM
   92 J=J+3
      K=K+1
      W(J)=XD(K)
      W(J+1)=YD(K)
      W(J+2)=WD(K)
      IF (K-M) 92,93,93
   93 I=W(J+2)+HLF
   94 IF (K-I) 95,95,96
   96 XD(K)=XD(MM)
      YD(K)=YD(MM)
      WD(K)=WD(MM)
      K=K-1
      MM=MM-1
      GO TO 94
   95 XD(K)=W(J)
      YD(K)=W(J+1)
      WD(K)=0.0D0
      K=K-1
      J=J-3
      IF (J) 91,91,93
   91 CALL VB06AD (M,N,XD,YD,WD,RD,XN,FN,GN,DN,THETA,IABS(IPRINT),W)
   90 RETURN
      END
      SUBROUTINE X04BAF(NOUT,REC)
C     MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986.
C
C     X04BAF writes the contents of REC to the unit defined by NOUT.
C
C     Trailing blanks are not output, except that if REC is entirely
C     blank, a single blank character is output.
C     If NOUT.lt.0, i.e. if NOUT is not a valid Fortran unit identifier,
C     then no output occurs.
C
C     .. Scalar Arguments ..
      INTEGER           NOUT
      CHARACTER*(*)     REC
C     .. Local Scalars ..
      INTEGER           I
C     .. Intrinsic Functions ..
      INTRINSIC         LEN
C     .. Executable Statements ..
      IF (NOUT.GE.0) THEN
C        Remove trailing blanks
         DO 20 I = LEN(REC), 2, -1
            IF (REC(I:I).NE.' ') GO TO 40
   20    CONTINUE
C        Write record to external file
   40    WRITE (NOUT,FMT=7786) REC(1:I)
      END IF
      RETURN
C
 7786 FORMAT (A)
      END
* X04FTEXT
*UPTODATE X04AAFTEXT
      SUBROUTINE X04AAF(I,NERR)
C     MARK 7 RELEASE. NAG COPYRIGHT 1978
C     MARK 7C REVISED IER-190 (MAY 1979)
C     IF I = 0, SETS NERR TO CURRENT ERROR MESSAGE UNIT NUMBER
C     (STORED IN NERR1).
C     IF I = 1, CHANGES CURRENT ERROR MESSAGE UNIT NUMBER TO
C     VALUE SPECIFIED BY NERR.
C
C     *** NOTE ***
C     THIS ROUTINE ASSUMES THAT THE VALUE OF NERR1 IS SAVED
C     BETWEEN CALLS.  IN SOME IMPLEMENTATIONS IT MAY BE
C     NECESSARY TO STORE NERR1 IN A LABELLED COMMON
C     BLOCK /AX04AA/ TO ACHIEVE THIS.
C     IN THE VAX IMPLEMENTATION, THE SAVE STATEMENT IS USED.
C
C     .. SCALAR ARGUMENTS ..
      INTEGER I, NERR
C     ..
C     .. LOCAL SCALARS ..
      INTEGER NERR1
C     ..
      SAVE NERR1
      DATA NERR1 /6/
      IF (I.EQ.0) NERR = NERR1
      IF (I.EQ.1) NERR1 = NERR
      RETURN
      END
** END OF X04AAFTEXT
*UPTODATE X04ABFTEXT
      SUBROUTINE X04ABF(I,NADV)
C     MARK 7 RELEASE. NAG COPYRIGHT 1978
C     MARK 7C REVISED IER-190 (MAY 1979)
C      IF I = 0, SETS NADV TO CURRENT ADVISORY MESSAGE UNIT NUMBER
C     (STORED IN NADV1).
C     IF I = 1, CHANGES CURRENT ADVISORY MESSAGE UNIT NUMBER TO
C     VALUE SPECIFIED BY NADV.
C
C     *** NOTE ***
C     THIS ROUTINE ASSUMES THAT THE VALUE OF NADV1 IS SAVED
C     BETWEEN CALLS.  IN SOME IMPLEMENTATIONS IT MAY BE
C     NECESSARY TO STORE NADV1 IN A LABELLED COMMON
C     BLOCK /AX04AB/ TO ACHIEVE THIS.
C     IN THE VAX IMPLEMENTATION, THE SAVE STATEMENT IS USED.
C
C     .. SCALAR ARGUMENTS ..
      INTEGER I, NADV
C     ..
C     .. LOCAL SCALARS ..
      INTEGER NADV1
C     ..
      SAVE NADV1
      DATA NADV1 /6/
      IF (I.EQ.0) NADV = NADV1
      IF (I.EQ.1) NADV1 = NADV
      RETURN
      END
** END OF X04ABFTEXT
*UPTODATE X04BZZTEXT
      SUBROUTINE X04BZZ(IFMTR)
C     MARK 8 RELEASE. NAG COPYRIGHT 1980.
C
C     **************** MACHINE-DEPENDENT SUBROUTINE ****************
C
C     X04BZZ SETS UP A FORMAT SPECIFICATION IN THE INTEGER ARRAY
C     IFMTR, FOR PRINTING A REAL VARIABLE TO FULL MACHINE PRECISION
C     IN THE CENTRE OF A FIELD 30 CHARACTERS WIDE.
C
C     SUPPOSE  S+1  SIGNIFICANT DIGITS ARE TO BE PRINTED.  THEN THE
C     FORMAT SPECIFICATION MUST HAVE THE FORM
C     .               MX,1PEW.S,NX
C     WHERE  M, W, S, N  ARE INTEGERS SUCH THAT
C     .     W = S + 7
C     .     M + N = 30 - W
C     .     M = N  IF  S  IS ODD, OR  M = N + 1  IF  S  IS EVEN.
C     ALL UNUSED ELEMENTS OF IFMTR MUST BE FILLED WITH BLANKS.
C
C     IMPLEMENTORS MUST COMPLETE THE ASSIGNMENTS TO ELEMENTS
C     1, 7, 8, 10, 11 AND 13 OF IFMTR.  NOTE THAT 2 CHARACTER
C     POSITIONS ARE PROVIDED FOR  W  AND  S .  IF ONLY 1 CHARACTER
C     IS NEEDED, THE FIRST CHARACTER POSITION SHOULD BE SET TO A
C     BLANK.
C
C     PHILIP E. GILL, ENID M. R. LONG, WALTER MURRAY AND
C     SUSAN M. PICKEN
C     D.N.A.C., NATIONAL PHYSICAL LABORATORY, ENGLAND.
C     JULY 1978
C
C     **************************************************************
C
C     .. ARRAY ARGUMENTS ..
      INTEGER IFMTR(20)
C     ..
C     .. LOCAL SCALARS ..
      INTEGER I1, I2, I3, I4, I5, I6, I7, I8, I9, IBLANK, ICOMMA,ID, IP,
     * IPOINT, IX, IZERO, J
C     ..
      DATA ICOMMA /1H,/
      DATA ID /1HD/
      DATA IP /1HP/
      DATA IPOINT /1H./
      DATA IX /1HX/
      DATA IZERO /1H0/
      DATA I1 /1H1/
      DATA I2 /1H2/
      DATA I3 /1H3/
      DATA I4 /1H4/
      DATA I5 /1H5/
      DATA I6 /1H6/
      DATA I7 /1H7/
      DATA I8 /1H8/
      DATA I9 /1H9/
      DATA IBLANK /1H /
      IFMTR(1) = I4
      IFMTR(2) = IX
      IFMTR(3) = ICOMMA
      IFMTR(4) = I1
      IFMTR(5) = IP
      IFMTR(6) = ID
      IFMTR(7) = I2
      IFMTR(8) = I2
      IFMTR(9) = IPOINT
      IFMTR(10) = I1
      IFMTR(11) = I5
      IFMTR(12) = ICOMMA
      IFMTR(13) = I4
      IFMTR(14) = IX
      DO 20 J=15,20
         IFMTR(J) = IBLANK
   20 CONTINUE
      RETURN
C
C     END OF SUBROUTINE X04BZZ
C
      END
** END OF X04BZZTEXT
        SUBROUTINE ZERO(A,N)
C       ====================
C
C ZERO N BYTES OF A
C
        INTEGER*1 A(N)
        DO 1 I=1,N
1       A(I)=0
        RETURN
        END
C*LLGRAPH**********************************************************************
C
C  PLOT AMPLITUDES AND PHASES ALONG EACH LAYER LINE ON PLOTTER.
C  ADAPTED FROM DAA LATLINE LATTICE-LINE FITTING PACKAGE
C
C    RAC  2NOV82   AUTO-SCALING FOR AMPLITUDES
C    RAC  11NOV00  CONVERT TO IMSUBS2000/PLOT_2K
C
C  NOBS    - NO. OF OBSERVATIONS OF AMP AND PHASE.
C  FOBS    - ARRAY OF OBSERVED AMPS.
C  PHIOBS  - ARRAY OF OBSERVED PHASES.
C  ZSTAR   - ARRAY OF OBSERVED ZSTARS.
C  ZMIN    - MINIMUM Z IN PLOT.
C  ZMAX    - MAXIMUM Z IN PLOT.
C  FMAX    - MAXIMUM AMP IN PLOT.
C  LLNO    - LAYER LINE NO. This and following just annotate plot.
C  SPLL    - LAYER LINE SPACING IN RECIPROCAL ANGSTROMS
C  NBES    - BESSEL ORDER
C  TITLE   - TITLE FOR PLOT
C  LAST    - .TRUE. IF THIS IS LAST PLOT (ELSE= .FALSE.)
C
C    FSCALE HEIGHT=70MM,PSCALE HEIGHT=50MM.
C    ZSCALE IS NOW 170MM
C
C
	SUBROUTINE LLGRAPH(FOBS,PHIOBS,ZSTAR,ZMIN,ZMAX,FMAX,LLNO,SPLL,
     1   NBES,NOBS,TITLE,LAST)
	DIMENSION FOBS(1),PHIOBS(1),ZSTAR(1),TITLE(1)
CTSH	DIMENSION LINE(20)
	character	line*80
	DATA INIT/0/
	LOGICAL LAST
C
C	ZMAG is width of boxes, FMAG is height of amp box,
C	PMAG is height of phase box, GAP is gap between boxes, all in mm
C	DELZ is interval for resolution scale in reciprocal Angstroms
	ZMAG=170.
	FMAG=70.
	PMAG=50.
	GAP=5.
	DELZ=0.01
C
	IF(INIT.EQ.1) GO TO 5
	CALL P2K_OUTFILE('llplot.ps',9)
	INIT=1
5	ZRANG=ZMAX-ZMIN
	ZERO=ZMAG/2.
C
C	Draw title and amplitude box, setting plotting area to +/-90mm
        CALL P2K_HOME
        CALL P2K_GRID(90.,90.,1.0)
	CALL P2K_ORIGIN(-90.0,-90.0,0)
        CALL P2K_MOVE(10.,-10.,0.)
	CALL P2K_FONT('Courier'//CHAR(0),3.0)
	CALL P2K_STRING(TITLE,40,0.0)
	CALL P2K_MOVE(0.0,0.0,0.)
	CALL P2K_DRAW(0.0,FMAG,0.0)
	CALL P2K_DRAW(ZMAG,FMAG,0.0)
	CALL P2K_DRAW(ZMAG,0.0,0.0)
	CALL P2K_DRAW(0.0,0.0,0.0)
	CALL P2K_MOVE(ZERO,0.0,0.0)
	CALL P2K_DRAW(ZERO,FMAG,0.0)
C	Write layer line details
	POSNX=ZMAG-40.
	POSNY=FMAG-5.
	CALL P2K_MOVE(POSNX,POSNY,0.0)
CTSH	ENCODE(20,103,LINE) LLNO
103	FORMAT('LAYER LINE    ',I4)
	write(line,'(i4)') llno
	CALL P2K_STRING(LINE,20,0.0)
        POSNY=POSNY-5.
        CALL P2K_MOVE(POSNX,POSNY,0.0)
CTSH        ENCODE(20,104,LINE) NBES
104	FORMAT('BESSEL ORDER  ',I4)
	write(line,'(i4)') nbes
        CALL P2K_STRING(LINE,20,0.0)
	POSNY=POSNY-5.
        CALL P2K_MOVE(POSNX,POSNY,0.0)
	write(line,'(f6.5)') spll
CTSH        ENCODE(20,105,LINE) SPLL
105     FORMAT('ZSTAR       ',F6.5)
        CALL P2K_STRING(LINE,20,0.0)
C	Put scale on radial axis
C	Number of tick marks
	IZ=ZRANG/DELZ
	IZ2=IZ/2
	ZA=-IZ2*DELZ
C	mm spacing of tick marks
	DELX=ZMAG*DELZ/ZRANG
	XPOS=ZERO-IZ2*DELX
C
	DO 25 J=1,IZ
	  CALL P2K_MOVE(XPOS,0.0,0.0)
	  CALL P2K_DRAW(XPOS,2.0,0.0)
	  XCH=XPOS-4.0
	  CALL P2K_MOVE(XCH,-5.0,0.0)
CTSH	  ENCODE(6,26,LINE) ZA
	  write(line(1:24),'(f6.2)') zpos
	  CALL P2K_STRING(LINE,6,0.0)
	  ZA=ZA+DELZ
	  XPOS=XPOS+DELX
25	CONTINUE
26	FORMAT(F6.2)
	POSNX=ZMAG-2.
	POSNY=-3.
	CALL P2K_FONT('Courier'//CHAR(0),2.0)
	CALL P2K_MOVE(POSNX,POSNY,0.0)
	CALL P2K_STRING('RECIPROCAL',10,0.0)
	CALL P2K_MOVE(POSNX,POSNY-3.0,0.0)
	CALL P2K_STRING('ANGSTROMS',9,0.0)
C
C	Put scale for amplitudes
	CALL P2K_FONT('Courier'//CHAR(0),3.0)
	SCALE=FMAG/(1.05*FMAX)
	IA=ALOG10(1.05*FMAX)
	B=10.0**IA
	IC=FMAX*1.05/B
	DO 200 J=1,IC
	  F=J*B
	  YPOS=F*SCALE
	  ZA=0.
	  ZB=ZMAG
	  CALL P2K_MOVE(ZA,YPOS,0.0)
	  ZD=ZA+2.0
	  CALL P2K_DRAW(ZD,YPOS,0.0)
	  ZD=ZB-2.0
	  CALL P2K_MOVE(ZB,YPOS,0.0)
	  CALL P2K_DRAW(ZD,YPOS,0.0)
	  XPOS=ZB
	  CALL P2K_MOVE(XPOS,YPOS,0.0)
CTSH	  ENCODE(7,201,LINE) F
	  write(line(1:28),'(f7.1)') f
	  CALL P2K_STRING(LINE,7,0.0)
200	CONTINUE
201	FORMAT(F7.1)
C
C	Plot observed amplitudes
C
	CALL P2K_FONT('Courier'//CHAR(0),1.0)
	DO 50 J=1,NOBS
	  IF(FOBS(J).EQ.-999.) GO TO 50
	  XP=ZERO+ZSTAR(J)*ZMAG/ZRANG
	  YP=FOBS(J)*SCALE
	  CALL P2K_MOVE(XP,YP,0.0)
	  CALL P2K_CSTRING('X',1,0.0)
50	CONTINUE
	CALL P2K_FONT('Courier'//CHAR(0),3.0)
C
C	Draw phase box
C
	PMAG2 = PMAG/360.
	YPOS=FMAG+GAP+PMAG/2.0
	CALL P2K_MOVE(0.0,0.0,0.0)
	CALL P2K_ORIGIN(0.0,YPOS,0.0)
	YPOS=PMAG/2.0
	CALL P2K_MOVE(0.0,-YPOS,0.0)
	CALL P2K_DRAW(ZMAG,-YPOS,0.0)
	CALL P2K_DRAW(ZMAG,YPOS,0.0)
	CALL P2K_DRAW(0.0,YPOS,0.0)
	CALL P2K_DRAW(0.0,-YPOS,0.0)
	CALL P2K_MOVE(ZERO,-YPOS,0.0)
	CALL P2K_DRAW(ZERO,YPOS,0.0)
C	Put scale for phases
	DO 620 J=1,7
	  YPOS=(PMAG/8.)*J-PMAG/2.
	  CALL P2K_MOVE(ZA,YPOS,0.0)
	  ZD=ZA+2.0
	  CALL P2K_DRAW(ZD,YPOS,0.0)
	  ZD=ZB-2.0
	  CALL P2K_MOVE(ZB,YPOS,0.0)
	  CALL P2K_DRAW(ZD,YPOS,0.0)
620	CONTINUE
	DO 630 J=1,5
	  IANG=-180+(J-1)*90
	  XPOS=ZB+1.0
	  YPOS=IANG*PMAG2
	  CALL P2K_MOVE(XPOS,YPOS,0.0)
CTSH	  ENCODE(4,631,LINE) IANG
	  write(line(1:16),'(i4)') iang
	  CALL P2K_STRING(LINE,4,0.0)
630	CONTINUE
631	FORMAT(I4)
C
C	Plot observed phases
C
	CALL P2K_FONT('Courier'//CHAR(0),1.0)
	DO 500 J=1,NOBS
	  IF(PHIOBS(J).EQ.-999.) GO TO 500
	  XP=ZERO+ZSTAR(J)*ZMAG/ZRANG
	  YP=PHIOBS(J)*PMAG2
	  CALL P2K_MOVE(XP,YP,0.0)
	  CALL P2K_CSTRING('X',1,0.0)
500	CONTINUE
C
	CALL P2K_PAGE
	RETURN
C
	END
C**************************************************************
CTSH this has been replaced by intrinsic FLOAT
CTSH	real function floati(i2)
CTSHC*** function to replace floati for compatibility between big
CTSHC*** and little endian machines
CTSH	integer*2	i2
CTSH	integer*4	i4
CTSH	i4 = i2
CTSH	floati = float(i4)
CTSH	return
CTSH	end
